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Introduction 



The aim of these lecture notes is to provide an introduction to methods and 
techniques used in the numerical solution of simple (non-relativistic) quantum- 
mechanical problems, with special emphasis on atomic and condensed-matter 
physics. The practical sessions are meant to be a sort of "computational lab- 
oratory" , introducing the basic ingredients used in the calculation of materials 
properties at a much larger scale. The latter is a very important field of today's 
computational physics, due to its technological interest and potential applica- 
tions. 

The codes provided during the course are little more than templates. Stu- 
dents are expected to analyze them, to run them under various conditions, to 
examine their behavior as a function of input data, and most important, to 
interpret their output from a physical point of view. The students will be asked 
to extend or modify those codes, by adding or modifying some functionalities. 

For further insight on the theory of Quantum Mechanics, many excellent 
textbooks are available (e.g. Griffiths, Schiff, or the ever-green Dirac and Lan- 
dau). For further insight on the properly computational aspects of this course, 
we refer to the specialized texts quotes in the Bibliography section, and in 
particular to the book of Thijssen. 

0.1 About Software 

This course assumes some basic knowledge of how to write and execute simple 
programs, and how to plot their results. All that is needed is a fortran or C 
compiler and some visualization software. The target machine is a PC running 
Linux, but other operating systems can be used as well (including Mac OS-X 
and Windows), as long as the mentioned software is installed and working, and 
if you know how to use it in oractise. 

0.1.1 Compilers 

In order to run a code written in any programming language, we must first 
translate it into machine language, i.e. a language that the computer can 
understand. The translation is done by an interpreter or by a compiler: the 
former translates and immediately executes each instruction, the latter takes 
the file, produces the so-called object code that together with other object codes 
and with libraries is finally assembled into an executable file. Python, Java (or at 



1 



an higher level, Matlab, Mathematica) are examples of "interpreted" language. 
Fortran, C, C++ are "compiled" languages. 

Our codes are written in Fortran 90. This is a sophisticated and com- 
plex language offering dynamical memory management, arrays operations (e.g. 
matrix- vector products), modular and object-based structure. Fortran 90 how- 
ever can be as efficient as Fortran 77 and maintains a wide compatibility with 
existing Fortran 77 codes. It is worth mentioning that the first applications 
of computers to physics go back to well before the birth of modern computer 
languages like C++, python, or even C: there is still a large number of codes 
and libraries written in Fortran 77 (or even Fortran 66!) widely used in physics. 

Fortran 90 (or even Fortran 77, in this respect) is not a well known language. 
There are however many available resources (see for instance the web page 
mentioned in the bibliography section) and the codes themselves are very simple 
and make little usage of advanced language features. In any case, there are no 
objections if a student prefers to use a more widespread language like C. A 
version of all codes in C is also available (no warranty about the quality of the 
C code in terms of elegance and good coding practice). 

In all cases, we need a C or Fortran 90 compiler. In PCs running Linux, the 
C compiler gcc is basically part of the operating system and is always present. 
Recent versions of gcc also include a Fortran compiler, called gf ortran. If this 
is absent, or it is not easy to install it, one can download the free and rather 
reliable compiler, g9^| It is possible to install on Mac OS-X and on Windows 
either gcc with gf ortran or g95. 



0.1.2 Visualization Tools 

Visualization of data produced by the codes (wave functions, charge densities, 
various other quantities) has a central role in the analysis and understanding 
of the results. Code gnuplot can be used to make two-dimensional or three- 
dimensional plots of data or of analytical expressions, gnuplot is open-source 
software, available for all operating systems and usually found pre-installed on 
Linux PCs. An introduction to gnuplot, with many links to more resources, 
can be found here: |http: //ww w, gnuplot . info/help .html, 

Another software that can be used is xmgracd^] This is also open-source 
and highly portable, has a graphical user interface and thus it is easier to use 
than gnuplot, whose syntax is not always easy to remember. 



0.1.3 Mathematical Libraries 

The usage of efficient mathematical libraries is crucial in "serious" calculations. 
Some of the codes use routines from the BLAS^] (Basic Linear Algebra Sub- 
programs) library and from LAPACKQ (Linear Algebra PACKage). The latter 

: http:/ /www. g95.org 

2 http : / / plasma-gate . weizmann. ac .il / Grace 
3 http:/ /www. netlib.org/blas 
' http://www.netlib.org/lapack 
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is an important and well-known library for all kinds of linear algebra oper- 
ations: solution of linear systems, eigenvalue problems, etc. LAPACK calls 
BLAS routines for all CPU- intensive calculations. The latter are available in 
highly optimized form for many different architectures. 

BLAS and LAPACK routines are written in Fortran 77. BLAS and LA- 
PACK are often available in many operating systems and can be linked directly 
by the compiler by adding -llapack -lblas. In the case of C compiler, it 
may be needed to add an underscore (_) in the calling program, as in: dsyev_, 
dgemm_. This is due to different C-Fortran conventions for the naming of "sym- 
bols" (i.e. compiled routines). Note that the C compiler may also need -lm to 
link general mathematical libraries (i.e. operations like the square root). 

0.1.4 Pitfalls in C-Fortran interlanguage calls 

In addition to the above-mentioned potential mismatches between C and For- 
tran naming conventions, there are a few more pitfalls one has to be aware of 
when Fortran 77 routines are called by C (or vice versa). 

• Fortran passes pointers to subroutines and functions; C passes values. In 
order to call a Fortran routine from C, all C variables appearing in the 
call must be either pointers or arrays. 

• Indices of vectors and arrays start from in C, from 1 in Fortran (unless 
differently specified in array declaration or allocation). 

• Matrices in C are stored in memory row- wise, that is: a[i] [j+1] follows 
a[i] [j] in memory. In Fortran, they are stored column- wise (the other 
way round!): a(i+l, j) follows a(i, j) in memory. 

An additional problem is that C does not provide run-time allocatable matrices 
like Fortran does, but only fixed-dimension matrices and arrays of pointers. 
The former are impractical, the latter are not usable as arguments to pass to 
Fortran. It would be possible, using either non-standard C syntax, or using 
C++ and the new command, to define dynamically allocated matrices similar 
to those used in Fortran. We have preferred for our simple C codes to "simulate" 
Fortran-style matrices (i.e. stored in memory column-wise) by mapping them 
onto one-dimensional C vectors. 

We remark that Fortran 90 has a more advanced way of passing arrays to 
subroutines using "array descriptors" . The codes used in this course however do 
not make use of this possibility but use the old-style Fortran 77 way of passing 
arrays via pointers. 
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Chapter 1 

One-dimensional Schrodinger 
equation 



In this chapter we will start from the harmonic oscillator to introduce a general 
numerical methodology to solve the one-dimensional, time-independent Schro- 
dinger equation. The analytical solution of the harmonic oscillator will be first 
derived and described. A specific integration algorithm (Numerov) will be used. 
The extension of the numerical methodology to other, more general types of 
potentials does not present any special difficulty. 

For a particle of mass m under a potential V(x), the one-dimensional, time- 
independent Schrodinger equation is given by: 

where i/}(x) is the wave function, in general complex, and h is the Planck con- 
stant h divided by 2ir. In the following we are focusing on the discrete spectrum: 



the set of isolated energy values for which Eq.(l.l) has normalizable solutions, 
localized in space. 



1.1 The harmonic oscillator 

The harmonic oscillator is a fundamental problem in classical dynamics as well 
as in quantum mechanics. It represents the simplest model system in which 
attractive forces are present and is an important paradigm for all kinds of vi- 
brational phenomena. For instance, the vibrations around equilibrium positions 
of a system of interacting particles may be described, via an appropriate coor- 
dinate transformation, in terms of independent harmonic oscillators known as 
normal vibrational modes. The same holds in quantum mechanics. The study 
of the quantum oscillator allows a deeper understanding of quantization and of 
its effects and of wave functions of bound states. 

In this chapter we will first remind the main results of the theory of the 
harmonic oscillator, then we will show how to set up a computer code that 
allows to numerically solve the Schrodinger equation for the harmonic oscillator. 
The resulting code can be easily modified and adapted to a different (not simply 
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quadratic) interaction potential. This will allow to study problems that, unlike 
the harmonic oscillator, do not have a simple analytical solution. 



1.1.1 Units 

The Schrodinger equation for a one-dimensional harmonic oscillator is, in usual 
notations: 

where K the force constant (the force on the mass being F = —Kx, proportional 
to the displacement x and directed towards the origin). Classically such an 
oscillator has a frequency (angular frequency) 

oj = J-. (1.3) 
y m 

It is convenient to work in adimensional units. These are the units that will 
be used by the codes presented at the end of this chapter. Let us introduce 
adimensional variables £, defined as 

/mK\ 1/4 /mw\ 1/2 



(using Eq.(1.3) for oj), and e, defined as 

E 



(1.5) 

By inserting these variables into the Schrodinger equation, we find 

which is written in adimensional units. 
1.1.2 Exact solution 

One can easily verify that for large £ (such that e can be neglected) the solutions 



of Eq.(1.6) must have an asymptotic behavior like 

<K0 ~ Ce ±e/2 (1.7) 

where n is any finite value. The + sign in the exponent must however be 
discarded: it would give raise to diverging, non-physical solutions (in which the 
particle would tend to leave the £ = point, instead of being attracted towards 
it by the elastic force) . It is thus convenient to extract the asymptotic behavior 
and assume 

m = H(Oe~Z 2 / 2 (1.8) 

where H (£) is a well-behaved function for large £ (i.e. the asymptotic behavior 
is determined by the second factor / 2 ). In particular, H(£) must not grow 
like e^ , or else we fall back into a undesirable non-physical solution. 
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Under the assumption of Eq.(1.8), Eq.(1.6) becomes an equation for H(^): 

H"(0 - 2£H'(0 + (2e - 1)H(0 = 0. (1.9) 

It is immediate to notice that £q = 1/2, -ffo(£) = 1 is the simplest solution. 
This is the ground state, i.e. the lowest-energy solution, as will soon be clear. 

In order to find all solutions, we expand H (£) into a series (in principle an 
infinite one): 



(1.10) 



n=0 



we derive the series to find H' and H" , plug the results into Eq. (1.9) and regroup 
terms with the same power of £. We find an equation 



J2 [(n + 2)(n + l)Ai+2 + (2e 

n=0 



2n - I) A n ] C = 



(1.11) 



that can be satisfied for any value of £ only if the coefficients of all the orders 
are zero: 

(n + 2){n + l)A n+2 + (2e - 2n - 1)A» = 0. (1.12) 



Thus, once ^4o an< 3 A\ are given, Eq.(1.12) allows to determine by recursion the 
solution under the form of a power series. 

Let us assume that the series contain an infinite number of terms. For large 
n, the coefficient of the series behave like 



A n +2 
An. 



2 

> 

n 



that is: A n ^ 



1 



(1.13) 



Remembering that exp(£ ) = E n s /^'i whose coefficient also behave as in 
Eq.(1.13), we see that recursion relation Eq.(1.12) between coefficients produces 
a function H(t) that grows like exp(£ 2 ), that is, produces unphysical diverging 
solutions. 

The only way to prevent this from happening is to have in Eq.(1.12) all 
coefficients beyond a given n vanish, so that the infinite series reduces to a 
finite-degree polynomial. This happens if and only if 
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n + 



where n is a non-negative integer. 

Allowed energies for the harmonic oscillator are thus quantized: 



(1.14) 



n H — I hui n = 0, 1, 2, 



(1.15) 



The corresponding polynomials H n (£) are known as Hermite polynomials. H n (t) 
is of degree n in £, has n nodes, is even [H n {— £) = H n {£)\ for even n, odd 
[H n (— = —H n (t)] for odd n. Since e - ^ I 2 is node-less and even, the complete 
wave function corresponding to the energy E n : 



MO = H n (0^ 
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(1.16) 
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Figure 1.1: Wave functions and probability density for the quantum harmonic 
oscillator. 



has n nodes and the same parity as n. The fact that all solutions of the 
Schrodinger equation are either odd or even functions is a consequence of the 
symmetry of the potential: V(—x) = V(x). 
The lowest-order Hermite polynomials are 

H (0 = 1, ffi(£) = 2£, H 2 {Z) = 4£ 2 -2, ff 3 (0 = $f - 12£. (1.17) 

A graph of the corresponding wave functions and probability density is shown 
in fig. [13} 

1.1.3 Comparison with classical probability density 

The probability density for wave functions ip n ( x ) of the harmonic oscillator 
have in general n + 1 peaks, whose height increases while approaching the 
corresponding classical inversion points (i.e. points where V{x) = E). 

These probability density can be compared to that of the classical harmonic 
oscillator, in which the mass moves according to x(t) = xq sin(wt). The proba- 
bility p(x)dx to find the mass between x and x + dx is proportional to the time 
needed to cross such a region, i.e. it is inversely proportional to the speed as a 
function of x: 

dx 

p(x)dx oc . (1-18) 
v(x) 

Since v(t) = xquj cos(wt) = u^Jxq — Xq sin 2 (u;t), we have 

p(x) oc Ji (1.19) 

Ixi — x 2 



This probability density has a minimum for x = 0, diverges at inversion points, 
is zero beyond inversion points. 

The quantum probability density for the ground state is completely different: 
has a maximum for x = 0, decreases for increasing x. At the classical inversion 
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point its value is still ~ 60% of the maximum value: the particle has a high 
probability to be in the classically forbidden region (for which V(x) > E). 

In the limit of large quantum numbers (i.e. large values of the index n), 
the quantum density tends however to look similar to the quantum one, but it 
still displays the oscillatory behavior in the allowed region, typical for quantum 
systems. 



1.2 Quantum mechanics and numerical codes: some 
observations 

1.2.1 Quantization 

A first aspect to be considered in the numerical solution of quantum problems 
is the presence of quantization of energy levels for bound states, such as for 



instance Eq.(1.15) for the harmonic oscillator. The acceptable energy values 



E n are not in general known a priori. Thus in the Schrodinger equation (1.1) 
the unknown is not just ip(x) but also E. For each allowed energy level, or 
eigenvalue, E n , there will be a corresponding wave function, or eigenf unction, 

1p n {x). 

What happens if we try to solve the Schrodinger equation for an energy E 
that does not correspond to an eigenvalue? In fact, a "solution" exists for any 
value of E. We have however seen while studying the harmonic oscillator that 
the quantization of energy originates from boundary conditions, requiring no 
unphysical divergence of the wave function in the forbidden regions. Thus, if 
E is not an eigenvalue, we will observe a divergence of ip(x). Numerical codes 
searching for allowed energies must be able to recognize when the energy is not 
correct and search for a better energy, until it coincides - within numerical or 
predetermined accuracy - with an eigenvalue. The first code presented at the 
end of this chapter implements such a strategy. 



1.2.2 A pitfall: pathological asymptotic behavior 

An important aspect of quantum mechanics is the existence of "negative" ki- 
netic energies: i.e., the wave function can be non zero (and thus the probability 
to find a particle can be finite) in regions for which V(x) > E, forbidden ac- 



cording to classical mechanics. Based on (1.1) and assuming the simple case in 



which V is (or can be considered) constant, this means 

d 2 i> 



<h , -k^(x) (1.20) 



where k 2 is a positive quantity. This in turns implies an exponential behavior, 



with both ip(x) ~ exp(kx) and ip(x) ~ exp(— kx) satisfying (1.20). As a rule 
only one of these two possibilities has a physical meaning: the one that gives 
raise to a wave function that decreases exponentially at large \x\. 

It is very easy to distinguish between the "good" and the "bad" solution for 
a human. Numerical codes however are less good for such task: by their very 
nature, they accept both solutions, as long as they fulfill the equations. If even 
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a tiny amount of the "bad" solution (due to numerical noise, for instance) is 
present, the integration algorithm will inexorably make it grow in the classically 
forbidden region. As the integration goes on, the "bad" solution will sooner or 
later dominate the "good" one and eventually produce crazy numbers (or crazy 
NaN's: Not a Number). Thus a nice-looking wave function in the classically 
allowed region, smoothly decaying in the classically forbidden region, may sud- 
denly start to diverge beyond some limit, unless some wise strategy is employed 
to prevent it. The second code presented at the end of this chapter implements 
such a strategy. 



1.3 Numerov's method 

Let us consider now the numerical solution of the (time-independent) Schrd- 
dinger equation in one dimension. The basic assumption is that the equation 
can be discretized, i.e. written on a suitable finite grid of points, and integrated, 
i.e. solved, the solution being also given on the grid of points. 

There are many big thick books on this subject, describing old and new 
methods, from the very simple to the very sophisticated, for all kinds of dif- 
ferential equations and all kinds of discretizations and integration algorithms. 
In the following, we will consider Numerov's method (named after Russian as- 
tronomer Boris Vasilyevich Numerov) as an example of a simple yet powerful 
and accurate algorithm. Numerov's method is useful to integrate second-order 
differential equations of the general form 

d 2 v 

= -g{x)y{x) + s{x) (1.21) 

where g(x) and s(x) are known functions. Initial conditions for second-order 
differential equations are typically given as 

y(xo) = 2/0, y'(x ) = y' . (1.22) 



The Schrddinger equation (1.1) has this form, with g(x) = j%[E — V(x)] and 



s(x) = 0. We will see in the next chapter that also the radial Schrddinger 
equations in three dimensions for systems having spherical symmetry belongs 
to such class. Another important equation falling into this category is Poisson's 
equation of electromagnetism, 

^ = -4^(*) (1-23) 

where p(x) is the charge density. In this case g(x) = and s(x) = —&irp(x). 

Let us consider a finite box containing the system: for instance, — x mSuX < 
x < x max , with x max large enough for our solutions to decay to negligibly small 
values. Let us divide our finite box into N small intervals of equal size, Ax 
wide. We call Xi the points of the grid so obtained, yj = y(x{) the values of the 
unknown function y(x) on grid points. In the same way we indicate by g% and 
Si the values of the (known) functions g{x) and s{x) in the same grid points. In 
order to obtain a discretized version of the differential equation (i.e. to obtain 
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an equation involving finite differences), we expand y(x) into a Taylor series 
around a point x n , up to fifth order: 



Vn-l = Vn - y'n^X + ^(Ax) 

+0[(Ax ' 

Vn+1 = 

+0[(Ax 



l y>»(Ax 



1 „,///// 
120 



y n + y' n Ax + ^(A*) 2 + |^'(Ax) 3 + ^"(A*) 4 + ^'"(A*) 5 



If we sum the two equations, we obtain: 

y n+ i + yn-i = 2y n + y'n{Ax) 2 + ^C(Ax) 4 + 0[(Ax)% 



(1.24) 
(1.25) 
(1.26) 



Eq.(1.21) tells us that 

Vn — dnVn > Sn — Zn- 

The quantity z n above is introduced to simplify the notations. The following 
relation holds: 

z n +i + z n -i = 2z n + 4'(Ax) 2 + 0[(Ax) 4 ] (1.27) 

(this is the simple formula for discretized second derivative, that can be obtained 
in a straightforward way by Taylor expansion up to third order) and thus 

//// // ^n+l 1 2Zn 

Zn ~ (Ax) 2 



Vn 



^ + 0[(Ax) 2 ]. 



(1.28) 



By inserting back these results into Eq.(jl.25|) one finds 

Vn+1 



= 2y n - yn-i + (-g n y n + s n )(Ax) 2 
+ j 2 -(-g n+1 y n+ i + s n+1 - g n -\y n -\ + s n _i + 2g n y n - 2s n )(Ax) 2 
+0[(Axf] 

(1-29) 

and finally the Numerov's formula 



Vn+1 



i + g n +i^rr- 



2y r , 



1 



5o (M! 



~ Vn-l 
, (A*) 2 



i , (Ax) 2 
1 +5„_i^ T 2^ 



(1.30) 



+(s n+ i + Ws n + Sn-xY-^f + 0[{Axf] 

that allows to obtain y n +\ starting from y n and y n -i, and recursively the func- 
tion in the entire box, as long as the value of the function is known in the first 
two points (note the difference with "traditional" initial conditions, Eq.(1.22), 
in which the value at one point and the derivative in the same point is speci- 
fied). It is of course possible to integrate both in the direction of positive x and 
in the direction of negative x. In the presence of inversion symmetry, it will be 
sufficient to integrate in just one direction. 

In our case — Schrodinger equation — the s n terms are absent. It is convenient 
to introduce an auxiliary array f n , defined as 



(Ax) 2 2m 
/ n = l + 5 n — — — , where g n = - V(x n )\. 

l£ n 

Within such assumption Numerov's formula can be written as 

(12 - Wf n )y n - fn-iVn-i 



Vn+1 



fn 



+1 



(1.31) 



(1.32) 
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1.3.1 Code: harmonicO 

Code harmonicO . f9C^] (or harmonicO . solves the Schrodinger equation for 
the quantum harmonic oscillator, using the Numerov's algorithm above de- 
scribed for integration, and searching eigenvalues using the "shooting method" . 



The code uses the adimensional units introduced in (1.4). 

The shooting method is quite similar to the bisection procedure for the 
search of the zero of a function. The code looks for a solution ijj n (x) with a 
pre-determined number n of nodes, at an energy E equal to the mid-point of the 
energy range [E mia , E max \, i.e. E = (E max + E min )/2. The energy range must 
contain the desired eigenvalue E n . The wave function is integrated starting 
from x = in the direction of positive x; tt the same time, the number of nodes 
(i.e. of changes of sign of the function) is counted. If the number of nodes is 
larger than n, E is too high; if the number of nodes is smaller than n, E is too 
low. We then choose the lower half-interval [-Emm, E max = E], or the upper half- 
interval [-Emm = E,E m3bX ], respectively, select a new trial eigenvalue E in the 
mid-point of the new interval, iterate the procedure. When the energy interval 
is smaller than a pre-determined threshold, we assume that convergence has 
been reached. 

For negative x the function is constructed using symmetry, since ■0 n (— x) = 
(— l) n ip n (x). This is of course possible only because V(—x) = V(x), otherwise 
integration would have been performed on the entire interval. The parity of the 
wave function determines the choice of the starting points for the recursion. For 
n odd, the two first points can be chosen as yo = an d an arbitrary finite value 
for y\. For n even, yo is arbitrary and finite, y\ is determined by Numerov's 



formula, Eq.(1.32), with f\ = /_i and y\=y-Y- 



(12 - lO/o)^ h 

yi = 777 • (i-33) 

The code prompts for some input data: 

• the limit x max for integration (typical values: 5 -j- 10); 

• the number N of grid points (typical values range from hundreds to a few 
thousand); note that the grid point index actually runs from to N, so 
that Ax = x max /iV; 

• the name of the file where output data is written; 

• the required number n of nodes (the code will stop if you give a negative 
number). 

Finally the code prompts for a trial energy. You should answer in order to 
search for an eigenvalue with n nodes. The code will start iterating on the 
energy, printing on standard output (i.e. at the terminal): iteration number, 
number of nodes found (on the positive x axis only) , the current energy eigen- 
value estimate. It is however possible to specify an energy (not necessarily an 



1 http://www.fisica.uniud.it/%7Egiannozz/Corsi/MQ/Software/F90/harmonicO.f90 



http:/ /www.fisica.uniud.it/%7Egiannozz/Corsi/MQ/Software/C/harmonic0.c 
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eigenvalue) to force the code to perform an integration at fixed energy and see 
the resulting wave function. It is useful for testing purposes and to better un- 
derstand how the eigenvalue search works (or doesn't work). Note that in this 
case the required number of nodes will not be honored; however the integra- 
tion will be different for odd or even number of nodes, because the parity of n 
determines how the first two grid points are chosen. 

The output file contains five columns: respectively, x, tp(x), | -^(rc) | 2 , p c \(x) 
and V(x). p c \(x) is the classical probability density (normalized to 1) of the 
harmonic oscillator, given in Eq.(1.19). All these quantities can be plotted as a 
function of x using any plotting program, such as gnuplot, shortly described in 
the introduction. Note that the code will prompt for a new value of the number 
of nodes after each calculation of the wave function: answer -1 to stop the code. 
If you perform more than one calculation, the output file will contain the result 
for all of them in sequence. Also note that the wave function are written for 
the entire box, from — x max to x max . 

It will become quickly evident that the code "sort of" works: the results 
look good in the region where the wave function is not vanishingly small, but 
invariably, the pathological behavior described in Sec. (1.2.2) sets up and wave 
functions diverge at large \x\. As a consequence, it is impossible to normalize 
the ifi(x). The code definitely needs to be improved. The proper way to deal 
with such difficulty is to find an inherently stable algorithm. 



1.3.2 Code: harmonicl 

Code harmonicl . f 9C^](or harmonicl . cQ is the improved version of harmonicO 

that does not suffer from the problem of divergence at large x. 

Two integrations are performed: a forward recursion, starting from x = 0, 
and a backward one, starting from x max . The eigenvalue is fixed by the condition 
that the two parts of the function match with continuous first derivative (as 
required for a physical wave function, if the potential is finite). The matching 
point is chosen in correspondence of the classical inversion point, x c \, he. where 
V( x ci) = E. Such point depends upon the trial energy E. For a function 
defined on a finite grid, the matching point is defined with an accuracy that is 
limited by the interval between grid points. In practice, one finds the index icl 
of the first grid point x c = icl Ax such that V(x c ) > E; the classical inversion 
point will be located somewhere between x c — Ax and x c . 

The outward integration is performed until grid point icl, yielding a func- 
tion iPl(x) defined in [0, x c ]; the number n of changes of sign is counted in the 
same way as in harmonicO. If n is not correct the energy is adjusted (lowered 
if n too high, raised if n too low) as in harmonicO. We note that it is not 
needed to look for changes of sign beyond x c : in fact we know a priori that in 
the classically forbidden region there cannot be any nodes (no oscillations, just 
decaying solution). 

If the number of nodes is the expected one, the code starts to integrate 
inward from the rightmost points. Note the statement y(mesh) = dx: its only 

http:/ /www. fisica.uniud.it/%7Egiannozz/Corsi/MQ/Software/F90/harmonicl.f90 
4 http:/ /www.fisica.uniud.it/%7Egiannozz/Corsi/MQ/Software/C/harmonicl.c 
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goal is to force solutions to be positive, since the solution at the left of the 
matching point is also positive. The value dx is arbitrary: the solution is 
anyway rescaled in order to be continuous at the matching point. The code 
stopst the same index icl corresponding to x c . We thus get a function ijjn(x) 
defined in [x c ,x max ]. 

In general, the two parts of the wave function have different values in x c : 
i/jl(,x c ) and tpji(x c ). We first of all re-scale i/jr(x) by a factor iPl(xc)/iPr(x c ), 
so that the two functions match continuously in x c . Then, the whole function 
ip(x) is renormalized in such a way that / \^(x)\ 2 dx = 1. 

Now comes the crucial point: the two parts of the function will have in 
general a discontinuity at the matching point ijj' R (x c ) — ip' L (x c ). This difference 
should be zero for a good solution, but it will not in practice, unless we are 
really close to the good energy E = E n . The sign of the difference allows us 
to understand whether E is too high or too low, and thus to apply again the 
bisection method to improve the estimate of the energy. 

In order to calculate the discontinuity with good accuracy, we write the 
Taylor expansions: 

vU = vl - y'i L ^ + b" L {^? + o[(Ax) 3 } 

y? +1 = y? + y?Ax + M R {Axf + 0[(Axf] [L ^> 

For clarity, in the above equation i indicates the index icl. We sum the two 
Taylor expansions and obtain, noting that y\ = yf = and that y'( L = y'{ R = 
y'( = —giyi as guaranteed by Numerov's method: 

yf-i + y?+i = zyi + (y\ R - y'i L )^ - gi yi(Ax) 2 + o[(Axf] (1.35) 

that is 

y ,n _ # _ + & -V- a^ftw + ol(Axf] (136) 



or else, by using the notations as in (1.31) 



yf - y't = Vtl+Vf+1 ^ m)Vl + 0[(Axf] (1.37) 

In this way the code calculated the discontinuity of the first derivative. If the 
sign of y^ R — y'^ is positive, the energy is too high (can you give an argument for 
this?) and thus we move to the lower half-interval; if negative, the energy is too 
low and we move to the upper half interval. As usual, convergence is declared 
when the size of the energy range has been reduced, by successive bisection, to 
less than a pre-determined tolerance threshold. 

During the procedure, the code prints on standard output a line for each 
iteration, containing: the iteration number; the number of nodes found (on the 
positive x axis only); if the number of nodes is the correct one, the discontinuity 
of the derivative y[ R — y[ L (zero if number of nodes not yet correct); the current 
estimate for the energy eigenvalue. At the end, the code writes the final wave 
function (this time, normalized to l!) to the output file. 
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1.3.3 Laboratory 



Here are a few hints for "numerical experiments" to be performed in the com- 
puter lab (or afterward), using both codes: 

• Calculate and plot eigenfunctions for various values of n. It may be 
useful to plot, together with eigenfunctions or eigenfunctions squared, the 
classical probability density, contained in the fourth column of the output 
file. It will clearly show the classical inversion points. With gnuplot, e.g.: 

plot "filename" u 1:3 w 1, "filename" u 1:4 w 1 

(u = using, 1:3 = plot column3 vs column 1, w 1 = with lines; the second 
"filename 1 ' can be replaced by ""). 

• Look at the wave functions obtained by specifying an energy value not 
corresponding to an eigenvalue. Notice the difference between the results 
of harmonicO and harmonicl in this case. 

• Look at what happens when the energy is close to but not exactly an 
eigenvalue. Again, compare the behavior of the two codes. 

• Examine the effects of the parameters xmax, mesh. For a given Ax, how 
large can be the number of nodes? 

• Verify how close you go to the exact results (notice that there is a con- 
vergence threshold on the energy in the code). What are the factors that 
affect the accuracy of the results? 



Possible code modifications and extensions: 



Modify the potential, keeping inversion symmetry. This will require very 
little changes to be done. You might for instance consider a "double-well" 
potential described by the form: 



V(x) 



+ 1 



e,6>Q. 



(1.38) 



• Modify the potential, breaking inversion symmetry. You might consider 
for instance the Morse potential: 



V(x) = D 



—2ax 



- 2e~ ax + 1 



(1.39) 



widely used to model the potential energy of a diatomic molecule. Which 
changes are needed in order to adapt the algorithm to cover this case? 
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Chapter 2 

Schrodinger equation for 
central potentials 



In this chapter we will extend the concepts and methods introduced in the pre- 
vious chapter for a one-dimensional problem to a specific and very important 
class of three-dimensional problems: a particle of mass m under a central po- 
tential V(r), i.e. depending only upon the distance r from a fixed center. The 
Schrodinger equation we are going to study in this chapter is thus 



ip(r) = Eip(r) 



(2.1) 



The problem of two interacting particles via a potential depending only upon 
their distance, V(\r\ — r^l), e.g. the Hydrogen atom, reduces to this case (see 
the Appendix), with m equal to the reduced mass of the two particles. 

The general solution proceeds via the separation of the Schrodinger equation 
into an angular and a radial part. In this chapter we will be consider the 
numerical solution of the radial Schrodinger equation. A non-uniform grid is 
introduced and the radial Schrodinger equation is transformed to an equation 
that can still be solved using Numerov's method. 



2.1 Variable separation 

Let us introduce a polar coordinate system (r, 6, (/>), where 6 is the polar angle, 
4> the azimuthal one, and the polar axis coincides with the z Cartesian axis. 
After some algebra, one finds the expression of the Laplacian operator: 



w2 1 d 

r 2 dr 



,d_ 
dr 



+ 



1 



8 



r 2 sin 8 86 



sin 



86 



+ 



1 



8 2 



r 2 sin 2 6 8(f) 2 



(2.2) 



It is convenient to introduce the operator L = L x + L y + L z , the square of the 
angular momentum vector operator, L = — ihr x V. Both L and L 2 act only 
on angular variables. In polar coordinates, the explicit representation of L 2 is 



L 



-ti 2 



1 8 

sk[6d6 



sin 



86 



+ 



8 2 



sin 2 6 dcj) 2 



(2.3) 
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The Hamiltonian can thus be written as 



H 



2m r 2 dr 



,d_ 

dr 



+ 



L 2 



2mr 2 



+ V{r). 



(2.4) 



The term L 2 /2mr 2 also appears in the analogous classical problem: the radial 
motion of a mass having classical angular momentum L c i can be described by an 
effective radial potential V(r) = V(r) + L 2 l /2mr 2 , where the second term (the 
"centrifugal potential" ) takes into account the effects of rotational motion. For 
high L c i the centrifugal potential "pushes" the equilibrium positions outwards. 

In the quantum case, both L 2 and one component of the angular momentum, 
for instance L z : 

L z = -ih~ (2.5) 

commute with the Hamiltonian, so L 2 and L z are conserved and H, L 2 , L z have 
a (complete) set of common eigenfunctions. We can thus use the eigenvalues of 
I? and L z to classify the states. Let us now proceed to the separation of radial 



and angular variables, as suggested by Eq.(2.4). Let us assume 



il){r,6,<t>) = R(r)Y(9, 



(2.6) 



After some algebra we find that the Schrodinger equation can be split into an 
angular and a radial equation. The solution of the angular equations are the 
spherical harmonics, known functions that are eigenstates of both L 2 and of 



L z Y f 



mfiYi 



i(l + i)h 2 Y t 



im \ 



(2.7) 



> and m = —I, ...,£ are integer numbers). 
The radial equation is 



h 2 i_ d 

2m r 2 dr 



,dR 



di 



+ 



V(r) + 



h 2 £(£ + 1) 
2mr 2 



R ni (r) = E nl R nS {r). (2i 



In general, the energy will depend upon I because the effective potential does; 
moreover, for a given £, we expect a quantization of the bound states (if existing) 
and we have indicated with n the corresponding index. 
Finally, the complete wave function will be 



Ipnimir, <t>) = Rne(r)Y £m (9, (j>) 



(2.9) 



The energy does not depend upon m. As already observed, m classifies the 
projection of the angular momentum on an arbitrarily chosen axis. Due to 
spherical symmetry of the problem, the energy cannot depend upon the orien- 
tation of the vector L, but only upon his modulus. An energy level E n i will 
then have a degeneracy 2t + 1 (or larger, if there are other observables that 
commute with the Hamiltonian and that we haven't considered). 
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2.1.1 Radial equation 

The probability to find the particle at a distance between r and r + dr from the 
center is given by the integration over angular variables: 



\ipnim(.r,9,(j))\ 2 rd9rsm9d(j)dr = \R ni \ 2 r 2 dr = \xni?dr 

where we have introduced the auxiliary function x( r ) as 

x{ r ) = rR(r) 

and exploited the normalization of the spherical harmonics: 

\Y em (9,(/))\ 2 d6 sm9d0 = l 



(2.10) 



(2.11) 



(2.12) 



(the integration is extended to all possible angles). As a consequence the nor- 
malization condition for \ is 



\Xnt(r)\ 2 dr = 1 



(2.13) 



The function |x( r )| 2 can thus be directly interpreted as a radial (probability) 
density. Let us re- write the radial equation for x{ r ) instead of R(r). Its is 
straightforward to find that Eq.(2.8) becomes 

h 2 d 2 x 



2m dr 2 



+ 



V(r) + 



h z l(l + 1) 

2mr 2 



E 



X {r) = 0. 



(2.14) 



We note that this equation is completely analogous to the Schrodinger equa- 



tion in one dimension, Eq.(l.l), for a particle under an effective potential: 

h 2 e(£ + 1) 



V(r) = V(r) + 



2mr 2 



(2.15) 



As already explained, the second term is the centrifugal potential. The same 
methods used to find the solution of Eq.(l.l), and in particular, Numerov's 
method, can be used to find the radial part of the eigenfunctions of the energy. 



2.2 Coulomb potential 

The most important and famous case is when V(r) is the Coulomb potential: 

where e = 1.6021 x 10~ 19 C is the electron charge, Z is the atomic number 
(number of protons in the nucleus), eo = 8.854187817 x 10~ 12 in MKSA units. 
Physicists tend to prefer the CGS system, in which the Coulomb potential is 
written as: 

V(r) = -Zq 2 /r. (2.17) 
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In the following we will use q e = e / (47reo) so as to fall back into the simpler 



2 2 

wing we win use C 

CGS form. 

It is often practical to work with atomic units (a.u.): units of length are 
expressed in Bohr radii (or simply, "bohr"), a$: 



h 2 

ao — ~ 

m e qi 



0.529177 A = 0.529177 x 10~ 10 m, (2.18) 



while energies are expressed in Rydberg (Ry): 

1 Ry = = 13.6058 eV. (2.19) 

when m e = 9.11 x 10~ 31 Kg is the electron mass, not the reduced mass of 
the electron and the nucleus. It is straightforward to verify that in such units, 
n = 1, m e = 1/2, q 2 e = 2. 

We may also take the Hartree (Ha) instead or Ry as unit of energy: 

4 

m a 

lHa = 2Ry = = 27.212 eV (2.20) 

thus obtaining another set on atomic units, in which h = l,m e = l,q e = 1- 
Beware! Never talk about "atomic units" without first specifying which ones. 
In the following, the first set ("Rydberg" units) will be occasionally used. 

We note first of all that for small r the centrifugal potential is the dominant 
term in the potential. The behavior of the solutions for r — > will then be 
determined by 

(2.2, 

yielding x(r) ~ r e+1 , G r x( r ) ~ r ~ t - The second possibility is not physical 
because x( r ) is n °t allowed to diverge. 

For large r instead we remark that bound states may be present only if 
E < 0: there will be a classical inversion point beyond which the kinetic energy 
becomes negative, the wave function decays exponentially, only some energies 
can yield valid solutions. The case E > corresponds instead to a problem of 
electron-nucleus scattering, with propagating solutions and a continuum energy 
spectrum. In this chapter, the latter case will not be treated. 

The asymptotic behavior of the solutions for large r — > oo will thus be 
determined by 

d 2 x 2m e 



d ,,- h ,-m-) (2.22) 

yielding %(r) ~ exp(±fcr), where k = \/—2m e E/h~. The + sign must be dis- 
carded as unphysical. It is thus sensible to assume for the solution a form 
like 

oo 

x (r) = /+ V fcr A nr n (2-23) 

n=0 

which guarantees in both cases, small and large r, a correct behavior, as long 
as the series does not diverge exponentially. 
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2.2.1 Energy levels 



The radial equation for the Coulomb potential can then be solved along the 
same lines as for the harmonic oscillator, Sect ]l.l| The expansion (2.23) is 



introduced into (2.14), a recursion formula for coefficients A n is derived, one 



finds that the series in general diverges like exp(2/cr) unless it is truncated to a 
finite number of terms, and this happens only for some particular values of E: 



Z 2 m e qj 
n 2 If? 



Ry 



(2.24) 



where n > £ + 1 is an integer known as main quantum number. For a given £ 
we find solutions for n = I + 1, £ + 2, . . .; or, for a given n, possible values for £ 
are £ = 0, 1, . . . , n — 1. 

Although the effective potential appearing in Eq.(2.14) depend upon £, and 



the angular part of the wave functions also strongly depend upon £, the energies 



(2.24) depend only upon n. We have thus a degeneracy on the energy levels 



with the same n and different £, in addition to the one due to the 2^ + 1 possible 



values of the quantum number m (as implied in (2.8) where m does not appear). 
The total degeneracy (not considering spin) for a given n is thus 



]T(2£+1) 



n 



(2.25) 



2.2.2 Radial wave functions 

Finally, the solution for the radial part of the wave functions is 



where 



2Z r / 2m e E n 



(2.27) 
n ao y Ti" 

and L? n t^{x) are Laguerre polynomials of degree n — £ — 1. The coefficient has 
been chosen in such a way that the following orthonormality relations hold: 

Xrd(r)Xn'l(.r)dr = 5 nn >. (2.28) 

The ground state has n = 1 and £ = 0: Is in "spectroscopic" notation (2p is 
n = 2, £ = 1, 3d is n = 3, £ = 2, 4/ is n = 4, £ = 3, and so on). For the Hydrogen 
atom (Z = 1) the ground state energy is — lRy and the binding energy of the 
electron is 1 Ry (apart from the small correction due to the difference between 
electron mass and reduced mass). The wave function of the ground state is a 
simple exponential. With the correct normalization: 

z 3/2 

Vioo(r, 9, 0) = -^^e~ Zr / a °. (2.29) 
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The dominating term close to the nucleus is the first term of the series, 
Xni{ r ) ~ r +1 . The larger £, the quicker the wave function tends to zero when 
approaching the nucleus. This reflects the fact that the function is "pushed 
away" by the centrifugal potential. Thus radial wave functions with large £ do 
not appreciably penetrate close to the nucleus. 

At large r the dominating term is x( r ) ~ r n exp(—Zr /nao). This means 
that, neglecting the other terms, |Xn^( r )| 2 has a maximum about r = n 2 ao/Z. 
This gives a rough estimate of the "size" of the wave function, which is mainly 
determined by n. 

In Eq.(2.26) the polynomial has n — I — 1 degree. This is also the number of 
nodes of the function. In particular, the eigenfunctions with I = have n — 1 
nodes; those with I = n — 1 are node-less. The form of the radial functions can 
be seen for instance on the Wolfram Research sit^] or explored via the Java 
applet at Davidson Colleg^] 



2.3 Code: hydrogen_radial 

The code hydrogen_radial . f 9C^]or hydrogen.radial . cj^Jsolves the radial equa- 
tion for a one-electron atom. It is based on harmonic 1, but solves a slightly 
different equation on a logarithmically spaced grid. Moreover it uses a more 
sophisticated approach to locate eigenvalues, based on a perturbative estimate 
of the needed correction. 

The code uses atomic (Rydberg) units, so lengths are in Bohr radii (ao = 1), 
energies in Ry, fr 2 /(2m e ) = 1, q 2 = 2. 



2.3.1 Logarithmic grid 



The straightforward numerical solution of Eq.(2.14) runs into the problem of 
the singularity of the potential at r = 0. One way to circumvent this difficulty 
is to work with a variable-step grid instead of a constant-step one, as done until 
now. Such grid becomes denser and denser as we approach the origin. "Serious" 
solutions of the radial Schrodinger in atoms, especially in heavy atoms, invari- 
ably involve such kind of grids, since wave functions close to the nucleus vary on 
a much smaller length scale than far from the nucleus. A detailed description of 
the scheme presented here can be found in chap. 6 of The Hartree-Fock method 
for atoms, C. Froese Fischer, Wiley, 1977. 

Let us introduce a new integration variable x and a constant-step grid in 
x, so as to be able to use Numerov's method without changes. We define a 
mapping between r and x via 

x = x(r). (2.30) 

The relation between the constant-step grid spacing Ax and the variable-step 
grid spacing is given by 

Ax = x'(r)Ar. (2.31) 
http:/ /library. wolfram.com/webMathematica/Physics/Hydrogen.jsp 



http:/ / webphysics.davidson.edu/physlet_resources/cise_qm/html/hydrogenic. html 



http:/ /www. Hsica.uniud.it/%7Egiannozz/Corsi/MQ/Software/F90/hydrogen_radial.f90 



http:/ /www.fisica.uniud.it/%7Egiannozz/Corsi/MQ/Software/C/hydrogen_radial.c 
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We make the specific choice 



x(r) = log 



Zr 
a 



(note that with this choice x is adimensional) yielding 

Ar 

r 



Ax 



(2.32) 



(2.33) 



The Ar/r ratio remains thus constant on the grid of r, called logarithmic grid, 
so defined. 



There is however a problem: by transforming Eq.(2.14) in the new vari 



able x, a term with first derivative appears, preventing the usage of Numerov's 
method (and of other integration methods as well). The problem can be cir- 
cumvented by transforming the unknown function as follows: 



y(x) 



1 



X(r(x)) 



(2.34) 



It is easy to verify that by transforming Eq.(2.14) so as to express it as a 



function of x and y, the terms containing first-order derivatives disappear, and 
by multiplying both sides of the equation by r 3//2 one finds 



dx 2 



+ 



2 -^r 2 (E-V(r))-(l+iy 



y(x) 







(2.35) 



where V(r) = —Zq 2 /r for the Coulomb potential. This equation no longer 
presents any singularity for r = 0, is in the form of Eq.(1.21), with 



9{x) 



2m e 
IF 



r{x) 2 (E- V(r(x))) - [£ + 1 



(2.36) 



and can be directly solved using the numerical integration formulae Eqs.(1.31) 



and Eq.(1.32) and an algorithm very similar to the one of Sect. 1.3.2 



Subroutine dojnesh defines at the beginning and once for all the values of 
r, y/r, r 2 for each grid point. The potential is also calculated once and for all 
in init.pot. The grid is calculated starting from a minimum value x = —8, 
corresponding to Zr m \ n ~ 3.4 x 10~ 3 Bohr radii. Note that the grid in r does 
not include r = 0: this would correspond to x = — oo. The known analytical 
behavior for r — > and r — > oo are used to start the outward and inward 
recurrences, respectively. 



2.3.2 Improving convergence with perturbation theory 

A few words are in order to explain this section of the code: 
i = icl 

ycusp = (y(i-l)*f (i-l)+f (i+l)*y(i+l)+10.d0*f (i)*y(i)) / 12. dO 

dfcusp = f (i) * (y(i) /ycusp - l.dO) 

! eigenvalue update using perturbation theory 

de = dfcusp/ddxl2 * ycusp*ycusp * dx 
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whose goal is to give an estimate, to first order in perturbation theory, of the 
difference be between the current estimate of the eigenvalue and its final value. 

Reminder: icl is the index corresponding to the classical inversion point. 
Integration is made with forward recursion up to this index, with backward 
recursion down to this index, icl is thus the index of the matching point 
between the two functions. The function at the right is rescaled so that the 
total function is continuous, but the first derivative dy/dx will be in general 
discontinuous, unless we have reached a good eigenvalue. 

In the section of the code shown above, y(icl) is the value given by Nu- 
merov's method using either icl-1 or icl+1 as central point; ycusp is the value 
predicted by the Numerov's method using icl as central point. The problem 
is that ycusp^y (icl) . 

What about if our function is the exact solution, but for a different problem? 
It is easy to find what the different problem could be: one in which a delta func- 
tion, vqS(x — x c ), is superimposed at x c =x(icl) to the potential. The presence 
of a delta function causes a discontinuity (a "cusp") in the first derivative, as 
can be demonstrated by a limit procedure, and the size of the discontinuity is 
related to the coefficient of the delta function. Once the latter is known, we 
can give an estimate, based on perturbation theory, of the difference between 
the current eigenvalue (for the "different" potential) and the eigenvalue for the 
"true" potential. 

One may wonder how to deal with a delta function in numerical integration. 
In practice, we assume the delta to have a value only in the interval Ax centered 
on y(icl). The algorithm used to estimate its value is quite sophisticated. Let 



us look again at Numerov's formula (1.32): note that the formula actually 



provides only the product y(icl)f (icl). Prom this we usually extract y(icl) 
since f(icl) is assumed to be known. Now we suppose that f (icl) has a 
different and unknown value f cusp, such that our function satisfies Numerov's 
formula also in point icl. The following must hold: 

fcusp*ycusp = f (icl) *y (icl) 

since this product is provided by Numerov's method (by integrating from icl-1 
to icl+l), and ycusp is that value that the function y must have in order to 
satisfy Numerov's formula also in icl. As a consequence, the value of df cusp 
calculated by the program is just f cusp-f (icl) , or 5f. 

The next step is to calculate the variation 5V of the potential V(r) appearing 



{2m e /n T )r 2 5V. Since f(x) = 1 + g{x)(Ax) 2 /12, we have Sg = 12/(Ax) 2 <5/, 



in Eq.(2.35) corresponding to 5f. By differentiating Eq.(2.36) one finds Sg(x 



and thus 



Ti 2 1 12 



8V= oTir^Sf- ( 2 -37) 

2m e r 2 (Ax) 2 J v 1 

First-order perturbation theory gives then the corresponding variation of the 
eigenvalue: 

3e = {ip\5V\ip) = f \y(x)\ 2 r(x) 2 5Vdx. (2.38) 
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Note that the change of integration variable from dr to dx adds a factor r to 
the integral: 



r ^ 



10 



f(r)dr= I f{r{x)) d ^f^dx= / f{r{x))r{x)dx. (2.39) 
^ dx 



We write the above integral as a finite sum over grid points, with a single 
non-zero contribution coming from the region of width Ax centered at point 
x c =x(icl). Finally: 

Se = \y(x c )\ 2 r(x c ) 2 5VAx = - ~ ™ \y(x c )\ 2 Ax8f (2.40) 

is the difference between the eigenvalue of the current potential (i.e. with a 
superimposed Delta function) and that of the true potential. This expression 
is used by the code to calculate the correction de to the eigenvalue. Since in 
the first step this estimate may have large errors, the line 

e = max(min(e+de,eup) ,elw) 

prevents the usage of a new energy estimate outside the bounds [elw, eip] . As 
the code proceeds towards convergence, the estimate becomes better and better 
and convergence is very fast in the final steps. 

2.3.3 Laboratory 

• Examine solutions as a function of n and £; verify the presence of acci- 
dental degeneracy. 

• Examine solutions as a function of the nuclear charge Z. 



• Compare the numerical solution with the exact solution, Eq.(2.29), for 
the Is case (or other cases if you know the analytic solution). 

• Slightly modify the potential as defined in subroutine init_pot, verify 
that the accidental degeneracy disappears. Some suggestions: V{r) = 
—Zq 2 /r 1+& where 5 is a small, positive or negative, number; or add an 
exponential damping (Yukawa) V{r) = — Zq 2 exp(—Qr)/r where Q is a 
number of the order of 0.05 a.u.. 

• Calculate the expectation values of r and of 1/r, compare them with the 
known analytical results. 

Possible code modifications and extensions: 

• Consider a different mapping: r(x) = ro(exp(x) — 1), that unlike the one 
we have considered, includes r = 0. Which changes must be done in order 
to adapt the code to this mapping? 
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Chapter 3 

Scattering from a potential 



Until now we have considered the discrete energy levels of simple, one-electron 
Hamiltonians, corresponding to bound, localized states. Unbound, delocalized 
states exist as well for any physical potential (with the exception of idealized 
models like the harmonic potential) at sufficiently high energies. These states 
are relevant in the description of elastic scattering from a potential, i.e. pro- 
cesses of diffusion of an incoming particle. Scattering is a really important 
subject in physics, since what many experiments measure is how a particle is 
deflected by another. The comparison of measurements with calculated results 
makes it possible to understand the form of the interaction potential between 
the particles. In the following a very short reminder of scattering theory is 
provided; then an application to a real problem (scattering of H atoms by rare 
gas atoms) is presented. This chapter is inspired to Ch.2 (pp. 14-29) of the book 
of Thijssen. 



3.1 Short reminder of the theory of scattering 

The elastic scattering of a particle by another is first mapped onto the equivalent 
problem of elastic scattering from a fixed center, using the same coordinate 
change as described in the Appendix. In the typical geometry, a free particle, 
described as a plane wave with wave vector along the z axis, is incident on the 
center and is scattered as a spherical wave at large values of r (distance from the 
center). A typical measured quantity is the differential cross section da(Q)/dQ, 
i.e. the probability that in the unit of time a particle crosses the surface in the 
surface element dS = r 2 d£l (where Q is the solid angle, d£l = sinOdOdcf), where 
6 is the polar angle and <p the azimuthal angle). Another useful quantity is the 
total cross section a to t = f(da(£l)/d£l)d£l. For a central potential, the system 
is symmetric around the z axis and thus the differential cross section does not 
depend upon <p. The cross section depends upon the energy of the incident 
particle. 

Let us consider a solution having the form: 

^(r) = e ikr + l^l e ikr (3.1) 
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with k = (0,0, k), parallel to the z axis. This solution represents an incident 
plane wave plus a diffused spherical wave. It can be shown that the cross section 
is related to the scattering amplitude f(9): 



dQ 



dQ = \f{9)\ 2 d£l = |/(6*)| 2 sin 



(3.2) 



Let us look for solutions of the form (3.1 ). The wave function is in general given 
by the following expression: 



oo I 



</>(r) = £ £ A 



lm 



Xl{r 



Yt 



(3.3) 



Z=0 m=-l 

which in our case, given the symmetry of the problem, can be simplified as 



OO / \ 

i=o r 



(3.4) 



The functions Xi( r ) are solutions for (positive) energy E = h 2 k 2 /2m with an- 
gular momentum I of the radial Schrodinger equation: 



2m dr 2 



+ 



E-V(r) 



h 2 l(i + i) 

2mr 2 



Xi{r) = 0. 



(3.5) 



It can be shown that the asymptotic behavior at large r of xi( r ) is 

Xi{f) — kr [ji(kr) cos Si — ni(kr) sin 5i] (3-6) 

where the ji and n\ functions are the well-known spherical Bessel functions, 
respectively regular and divergent at the origin. These are the solutions of the 
radial equation for Ri(r) = Xi{ r )/ r i n absence of the potential. The quantities 
Si are known as phase shifts. We then look for coefficients of Eq.(3.4) that yield 
the desired asymptotic behavior (3.1). It can be demonstrated that the cross 
section can be written in terms of the phase shifts. In particular, the differential 
cross section can be written as 



da _ 1 

dtt ~ ¥ 



£(2/ + l)e i&l sin S t Pi (cos I 

1=0 



(3.7) 



while the total cross section is given by 



47T 



= ^2 £(2/ + l)sin 2 ^ 



1=0 



'2mE 
1^' 



(3.8) 



Note that the phase shifts depend upon the energy. The phase shifts thus con- 
tain all the information needed to know the scattering properties of a potential. 
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3.2 Scattering of H atoms from rare gases 



The total cross section crtot(E) for the scattering 
of H atoms by rare gas atoms was measured by 
Toennies et al, J. Chem. Phys. 71, 614 (1979). 
At the right, the cross section for the H-Kr sys- f 
tern as a function of the energy of the center of § 
mass. One can notice "spikes" in the cross sec- f 
tion, known as "resonances". One can connect 
the different resonances to specific values of the 
angular momentum /. 

The H-Kr interaction potential can be modelled quite accurately as a Lennard- 
Jones (LJ) potential: 




Energy [meVJ 



V(r) 



12 



a 



(3.9) 



where e = 5.9meV, a = 3.57A. The LJ potential is much used in molecular and 
solid-state physics to model interatomic interaction forces. The attractive r -6 
term describes weak van der Waals (or "dispersive" , in chemical parlance) forces 
due to (dipole-induced dipole) interactions. The repulsive r~ 12 term models the 
repulsion between closed shells. While usually dominated by direct electrostatic 
interactions, the ubiquitous van der Waals forces are the dominant term for the 
interactions between closed-shell atoms and molecules. These play an important 
role in molecular crystals and in macro- molecules. The LJ potential is the first 
realistic interatomic potential for which a molecular- dynamics simulation was 
performed (Rahman, 1965, liquid Ar). 



It is straightforward to find that the LJ potential as written in (3.9) has a 
minimum V m i n = —e for r = a, is zero for r = cr/2 1 / 6 = 0.89<7 and becomes 
strongly positive (i.e. repulsive) at small r. 



3.2.1 Derivation of Van der Waals interaction 

The Van der Waals attractive interaction can be described in semiclassical terms 
as a dipole-induced dipole interaction, where the dipole is produced by a charge 
fluctuation. A more quantitative and satisfying description requires a quantum- 
mechanical approach. Let us consider the simplest case: two nuclei, located 
in Ryi and R^, and two electrons described by coordinates ri and r 2 . The 
Hamiltonian for the system can be written as 



H 



2m 



Ti 2 



|ri — Ha\ 2m 

a 2 

+ 



|r 2 - Rb| 



(3.10) 



+ 



|n-Rfl| |t2-Ra| |ri-r 2 | |R-a - Kb\ 

where Vj indicates derivation with respect to variable r^, i = 1,2. Even this 
"simple" hamiltonian is a really complex object, whose general solution will be 
the subject of several chapters of these notes. We shall however concentrate on 
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the limit case of two Hydrogen atoms at a large distance R, with R = R/i — Rg. 
Let us introduce the variables xi = ri— R^, x 2 = r 2 — R_b- In terms of these new 
variables, we have H = Ha + Hb + AH(R), where Ha (Hb) is the Hamiltonian 
for a Hydrogen atom located in R^ (R_b), and AH has the role of "perturbing 
potential" : 

AH = - t q2 " l - l q2 " l + ] ^ + 4- (3-11) 

|xi + R| |x 2 -R| |xi-x 2 + R| R v ; 

Let us expand the perturbation in powers of 1/R. The following expansion, 
valid for R — > oo, is useful: 

1 1 R x 3(R • x) 2 - x 2 R 2 

+ ~ • (3-12) 



|x + R| R R 3 R 5 
Using such expansion, it can be shown that the lowest nonzero term is 

A ff .f(x 1 . X2 -3<^f^>). (3.13) 

The problem can now be solved using perturbation theory. The unperturbed 
ground-state wavefunction can be written as a product of Is states centered 
around each nucleus: $o( x i> x 2) = V , is( x i)V'is( x 2) (it should actually be anti- 
symmetrized but in the limit of large separation it makes no difference.). It is 
straightforward to show that the first-order correction, A^E = (<3?o| AH\$q), 
to the energy, vanishes because the ground state is even with respect to both xi 
and X2. The first nonzero contribution to the energy comes thus from second- 
order perturbation theory: 

4 >o Ei ~ E ° 

where are excited states and Ei the corresponding energies for the unper- 
turbed system. Since AH oc R~ 3 , it follows that A^E = -C 6 /R 6 , the well- 
known behavior of the Van der Waals interaction^ The value of the so-called 



Cq coefficient can be shown, by inspection of Eq.(3.14), to be related to the 
product of the polarizabilities of the two atoms. 



3.3 Code: crossection 

Code crossection. f9(j^J or crossection. calculates the total cross section 
atot(E) and its decomposition into contributions for the various values of the 
angular momentum for a scattering problem as the one described before. 

The code is composed of different parts. It is always a good habit to verify 
separately the correct behavior of each piece before assembling them into the 
final code (this is how professional software is tested, by the way). The various 
parts are: 

1 Note however that at very large distances a correct electrodynamical treatment yields a 

diff erent behavior: AE oc —R~ 7 . 

http:/ /www. fisica.uniud.it/%7Egiannozz/Corsi/MQ/Software/F90/crossection.f90 
http:/ /www.fisica.uniud.it/%7Egiannozz/Corsi/MQ/Software/C/crossection.c 
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1. Solution of the radial Schrodinger equation, Eq.(3.5), with the Lennard- 



Jones potential (3.9) for scattering states (i.e. with positive energy). One 
can simply use Numerov's method with a uniform grid and outwards 
integration only (there is no danger of numerical instabilities, since the 
solution is oscillating). One has however to be careful and to avoid the 
divergence (as ~ r~ 12 ) for r — > 0. A simple way to avoid running into 
trouble is to start the integration not from r m i n = but from a small but 
nonzero value. A suitable value might be r m j ra ~ 0.5a, where the wave 
function is very small but not too close to zero. The first two points can be 
calculated in the following way, by assuming the asymptotic (vanishing) 
form for r — > 0: 



„ 2mea 12 . I \2mea 12 . 

X (r) ^ ^2-712 *( r ) =*■ X(H — exp I — W 2 r | (3.15) 



(note that this procedure is not very accurate because it introduces and 
error of lower order, i.e. worse, than that of the Numerov's algorithm. 
In fact by assuming such form for the first two steps of the recursion, 
we use a solution that is neither analytically exact, nor consistent with 
Numerov's algorithm). 

The choice of the units in this code is (once again) different from that of 
the previous codes. It is convenient to choose units in which h 2 /2m is a 
number of the order of 1. Two possible choices are meV/A 2 , or meV/cr 2 . 
This code uses the former choice. Note that m here is not the electron 
mass! it is the reduced mass of the H-Kr system. As a first approximation, 
m is here the mass of H. 

2. Calculation of the phase shifts 61(E). Phase shifts can be calculated by 
comparing the calculated wave functions with the asymptotic solution at 
two different values r\ and r 2 , both larger than the distance r max beyond 
which the potential can be considered to be negligible. Let us write 

Xi(ri) = Akn \ji(kn) cos Si - ni(kri) sin Si] (3.16) 
Xi(r 2 ) = Akr 2 [ji(kr 2 ) cos 5i - ni(kr 2 ) sin 5i], (3.17) 

from which, by dividing the two relations, we obtain an auxiliary quantity 
K 

R = r 2 Xi{ri) = jijkn) - ni(kri)\,anSi 
nxifa) ji(kr 2 ) - ni(kr 2 ) tan 5i 
from which we can deduce the phase shift: 

t*nS l = ?J l ^ r2 \- jl ^ X \ . (3.19) 
Kni(kr 2 ) - ni{kr\) 

The choice of n and r 2 is not very critical. A good choice for r\ is about 
at r max . Since the LJ potential decays quite quickly, a good choice is 
T\ = r max = 5a. For r 2 , the choice is r 2 = r\ + A/2, where A = 1/k is the 
wave length of the scattered wave function. 
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3. Calculation of the spherical Bessel functions ji and n\. The analytical 
forms of these functions are known, but they become quickly unwieldy for 
high I. One can use recurrence relation. In the code the following simple 
recurrence is used: 

21 + 1 

z t+ i{x) = z t (x) - z = j,n (3.20) 

x 

with the initial conditions 

cosx . , , sin a; , > sinx . . cosx 



(x) = , jo(x) = ; n-i(x) = , n (x) 



XX XX 

(3-21) 

Note that this recurrence is unstable for large values of I: in addition to 
the good solution, it has a "bad" divergent solution. This is not a serious 
problem in practice: the above recurrence should be sufficiently stable up 
to at least I = 20 or so, but you may want to verify this. 



4. Sum the phase shifts as in Eq.(3.8) to obtain the total cross section and 
a graph of it as a function of the incident energy. The relevant range of 
energies is of the order of a few meV, something like 0.1 < E < 3H-4 
meV. If you go too close to E = 0, you will run into numerical difficulties 
(the wave length of the wave function diverges) . The angular momentum 
ranges from to a value of l max to be determined. 

The code writes a file containing the total cross section and each angular mo- 
mentum contribution as a function of the energy (beware: up to l ma x = 13; for 
larger I lines will be wrapped). 

3.3.1 Laboratory 

• Verify the effect of all the parameters of the calculation: grid step for in- 
tegration, r m i n , n = r max , r2, lmax- In particular the latter is important: 
you should start to find good results when l max > 6. 

• Print, for some selecvted values of the energy, the wave function and its 
asymptotic behavior. You should verify that they match beyond r max . 

• Observe the contributions of the various I and the peaks for I = 4, 5, 6 
(resonances). Make a plot of the effective potential as a function of I: can 
you see why there are resonances only for a few values of /? 

Possible code modifications and extensions: 

• Modify the code to calculate the cross section from a different potential, 
for instance, the following one: 

V(r) = -Aexp -(r-r ) 2 , r < r max ; V(r > r max ) = 

What changes do you need to apply to the algorithm, in addition to 
changing the form of the potential? 
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• In the limit of a weak potential (such as e.g. the potential introduced 
above) , the phase shifts 5i are well approximated by the Born approxima- 
tion: 

Si ~ K-k / r 2 jf(kr)V(r)dr, 

h Jo 

(k = ^2mE/K). Write a simple routine to calculate this integral, compare 
with numerical results. 
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Chapter 4 

The Variational Method 



The exact analytical solution of the Schrodinger equation is possible only in 
a few cases. Even the direct numerical solution by integration is often not 
feasible in practice, especially in systems with more than one particle. There are 
however extremely useful approximated methods that can in many cases reduce 
the complete problem to a much simpler one. In the following we will consider 
the variational principle and its consequences. This constitutes, together with 
suitable approximations for the electron-electron interactions, the basis for most 
practical approaches to the solution of the Schrodinger equation in condensed- 
matter physics. 

4.1 Variational Principle 

Let us consider a Hamiltonian H and a function ip, that can be varied at will 
with the sole condition that it stays normalized. One can calculate the expec- 
tation value of the energy for such function (in general, not an eigenfunction of 



where v represents all the integration coordinates. 

The variational principle states that functions ip for which (H) is stationary — 
i.e. does not vary to first order in small variations of tp — are the eigenfunctions 
of the energy. In other words, the Schrodinger equation is equivalent to a sta- 
tionarity condition. 

4.1.1 Demonstration of the variational principle 

Since an arbitrary variation 5ip of a wave function in general destroys its nor- 
malization, it is convenient to use a more general definition of expectation value, 
valid also for non-normalized functions: 



H): 




(4.1) 



(H) = 



$i)*H4>dv 



(4.2) 



/ ip*ipdv 
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By modifying ip as ip + Sip, the expectation value becomes 

jxr + sr)H(ip + 5ip)dv 

{ ' [ ' J(iP* + SiP*){iP + SiP) dv 

J^*HiPdv + J5ip*Hipdv + Jip*H5^dv 
J ip*ip dv + / 5ip*ip dv + / ip*Sip dv 

J iP*Hipdv + J 5ip*Hipdv + J ip*H5ipdv^j x 

i / [sripdv jr^dv \ 

J ip*ipdv \ J ip*ipdv ftp*ipdv J 

where second-order terms in dip have been omitted and we have used the ap- 
proximation 1/(1 + 2;) ~ 1 — x, valid for x « 1. By omitting again higher-order 
terms: 

s _ J Sip* Hip dv + Jip*H6ip*dv /J5ip*ipdv + Jip*5ipdv\ 

J ip*ipdv J ip*ipdv V / ip*ipdv j tp*ipdv J 

One of the two terms in parentheses is the complex conjugate of the other; the 
same holds for the two remaining terms, because H is a hermitian operator, 
satisfying 

' a*Hbdv = ( f b*Hadv) (4.5) 



for any pair of functions a and b. We can thus simplify the above expression as 

r)m fj5iP*Hipdv \ /rn fj5iP*ipdv \ 

5{H) = J „ , , , h c.c. - (H)[ , r + c.c. . (4.6) 

x 1 \ J ip*ip dv J \J 4>*ip dv J v 1 

Let us now assume that ip is such that (H) is stationary with respect to any 
variation of it: then 5(H) = 0, i.e. 



Sip* [H - (H)]ipdv + c.c. = (4.7) 

for an arbitrary variation Sip. This implies 

[H - (H)]iP = (4.8) 
that is, ip is a solution of the Schrodinger equation: 

Hip = Eip (4.9) 



4.1.2 Alternative demonstration of the variational principle 

A different and more general way to demonstrate the variational principle, 
which will be useful later, is based upon Lagrange multipliers method. This 
method deals with the problem of finding stationarity conditions for an integral 
Iq while keeping at the same time constant other integrals I\ . . .Ik- One can 
solve the equivalent problem 




(4.10) 
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where A& are constants to be determined (Lagrange multipliers). In our case we 
have 

J = J ip*H^jdv (4.11) 
h = Jtf>*tf>dv (4.12) 



and thus we assume 



S(I + Xh) = (4.13) 



where A must be determined. By proceeding like in the previous section, we 
find 

5I = J 6ip*Hif)dv + c.c. (4.14) 
<5Ji = J Sip*ipdv + c.c. (4.15) 
and thus the condition to be satisfied is 

<5(J + AJi) = J Sip* [H + A]V> dv + c.c. = (4.16) 

that is 

Hip = -Xtf; (4.17) 

i.e. the Lagrange multiplier is equal, apart from the sign, to the energy eigen- 
value. Again we see that states whose expectation energy is stationary with 
respect to any variation in the wave function are the solutions of the Schrodinger 
equation. 

4.1.3 Ground state energy 

Let us consider the eigenfunctions ip n of a Hamiltonian H, with associated 
eigenvalues (energies) E n : 

H^ n = E n ^ n . (4.18) 

Let us label the ground state with n = and the ground-state energy as Eq. 
Let us demonstrate that for any different function ip, we necessarily have 

{H) ~ fWdv ~ E °- (419) 

In order to demonstrate it, let us expand ip using as basis energy eigenfunctions 
(this is always possible because energy eigenfunctions are a complete orthonor- 
mal set): 

Y> = 5>„Vn (4.20) 



Then one finds 



itt\ J2n\ c n\ 2 E n J2n\ c n\ 2 {E n Eq) fAn-i\ 

\ H I = v \c |2 = + V \r I 2 ( ^ 



34 



This demonstrates Eq.(4.19), since the second term is either positive or zero, 
as E n > Eq by definition of ground state, 

This simple result is extremely important: it tells us that any function -0, 
yields for the expectation energy an upper estimate of the energy of the ground 
state. If the ground state is unknown, an approximation to the ground state 
can be found by varying ijj inside a given set of functions and looking for the 
function that minimizes (H). This is the essence of the variational method. 

4.1.4 Variational method in practice 

One identifies a set of trial wave functions ip(v;ai, . . . ,a r ), where v are the 
variables of the problem (coordinates etc), aj, i = 1, . . . , r are parameters. The 
energy eigenvalue will be a function of the parameters: 

E( ai ,...,a r ) = JiP*H^dv (4.22) 

The variational method consists in looking for the minimum of E with respect 
to a variation of the parameters, that is, by imposing 

dE dE 
dai da r 

The function i/j satisfying these conditions with the lowest E is the function 
that better approximates the ground state, among the considered set of trial 
functions. 

It is clear that a suitable choice of the trial functions plays a crucial role 
and must be carefully done. 

4.2 Secular problem 

The variational method can be reduced to an algebraic problem by expanding 
the wave function into a finite basis of functions, and applying the variational 



principle to find the optimal coefficients of the expansion. Based on Eq. (4.10), 
this means calculating the functional (i.e. a "function" of a function): 

ijj*H^dv-e f Tp*ijjdv (4.24) 



and imposing the stationary condition on G[ip]. Such procedure produces an 
equation for the expansion coefficients that we are going to determine. 

It is important to notice that our basis is formed by a finite number N of 
functions, and thus cannot be a complete system: in general, it is not possible 
to write any function tp (including exact solutions of the Schrodinger equation) 
as a linear combination of the functions in this basis set. What we are going to 
do is to find the ip function that better approaches the true ground state, among 
all functions that can be expressed as linear combinations of the iV chosen basis 
functions. 
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4.2.1 Expansion into a basis set of orthonormal functions 

Let us assume to have a basis of N functions bi, between which orthonormality 
relations hold: 

{bi\bj) = [ b*bjdv = Sij (4.25) 



Let us expand the generic ip in such basis: 

N 



1p = Yl ( 4 - 26 ) 



i=l 



By replacing Eq.(4.26) into Eq.(4.24) one can immediately notice that the latter 
takes the form 



G(ci,...,c N ) = ^2c*CjHij - e^2c*Cj5, 



= X''>;i", "V' (4-27) 

ij 

where we have introduced the matrix elements Hij: 

Hij = (bi\H\bj) = J b*Hbjdv (4.28) 

Since both H and the basis functions are given, Hij is a known square matrix of 
numbers. The hermiticity of the Hamiltonian operator implies that such matrix 
is hermitian: 

Hji = H* 3 (4.29) 

(i.e. symmetric if all elements are real). According to the variational method, 
let us minimize Eq. (4.27) with respect to the coefficients: 

sh° < 4 - 30 > 

This produces the condition 

J2(Hij - eS^Cj = (4.31) 
j 

If the derivative with respect to complex quantities bothers you: write the 
complex coefficients as a real and an imaginary part c& = Xk + iyk, require that 
derivatives with respect to both xt and yk are zero. By exploiting hermiticity 
you will find a system 

w k + w* k = 

-iW k + iW£ = 

where Wf, = J2j(^kj — £&kj) c ji that allows as only solution W k = 0. 

We note that, if the basis were a complete (and thus infinite) system, this 
would be the form of the Schrodinger equation. We have finally demonstrated 
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that the same equations, for a finite basis set, yield the best approximation to 
the true solution according to the variational principle. 



Eq.(4.31) is a system of N algebraic linear equations, homogeneous (there 
are no constant term) in the N unknown Cj. In general, this system has only 
the trivial solution Cj = for all coefficients, obviously unphysical. A non- 
zero solution exists if and only if the following condition on the determinant is 
fulfilled: 

det|fl*y -e5ij\ =0 ( 4 - 32 ) 

Such condition implies that one of the equations is a linear combination of the 
others and the system has in reality N — 1 equations and N unknowns, thus 
admitting non-zero solutions. 



Eq.(4.32) is known as secular equation. It is an algebraic equation of de- 
gree N in e (as it is evident from the definition of the determinant, with the 
main diagonal generating a term e N , all other diagonals generating lower-order 



terms), that admits ./V roots, or eigenvalues. Eq.(4.31) can also be written in 
matrix form 

He = ec (4.33) 

where H is here the NxN matrix whose matrix elements are Hij , c is the vector 
formed with a components. The solutions c are also called eigenvectors. For 
each eigenvalue there will be a corresponding eigenvector (known within a mul- 
tiplicative constant, fixed by the normalization). We have thus TV eigenvectors 
and we can write that there are N solutions: 

rpk = J2 C ikbi, k = l,...,N (4.34) 

i 

where Cik is a matrix formed by the N eigenvectors (written as columns and 
disposed side by side): 

Hifj k = ekipk, (4.35) 
that is, in matrix form, taking the i— th component, 

(fty fc )i = E H H C * = tkC ik . (4.36) 

3 



Eq.(4.33) is a common equation in linear algebra and there are standard 
methods to solve it. Given a matrix H, it is possible to obtain, using standard 
library routines, the C matrix and a vector e of eigenvalues. 

The solution process is usually known as diagonalization. This name comes 



from the following important property of C. Eq.(4.34) can be seen as a trans- 
formation of the ./V starting functions into another set of N functions, via a 
transformation matrix. It is possible to show that if the hi functions are or- 
thonormal, the ipk functions are orthonormal as well. Then the transformation 
is unitary, i.e. 

^C*jC ik = S jk (4.37) 



or, in matrix notations, 



(C-%- = C; = 4 (4.38) 
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that is, the inverse matrix is equal to the conjugate of the transpose matrix, 
known as adjoint matrix. The matrix C having such property is also called a 
unitary matrix. 

Let us consider now the matrix product C~ l HC and let us calculate its 
elements: 

(C l HC)kn = 1 )kiHijCj n 
ij 

En* \ u r 1 

i j 

= ^ C* k e n Ci n 

i 

= £n ^ C* k Ci n 

i 

= e n 5 kn (4.39) 

where the preceding results have been used. The transformation C reduces H 
to a diagonal matrix, whose non-zero N elements are the eigenvalues. We can 
thus see our eigenvalue problem as the search for a transformation that brings 
from the original basis to a new basis in which the H operator has a diagonal 
form, that is, it acts on the elements of the basis by simply multiplying them 
by a constant (as in Schrodinger equation). 



4.3 Plane-wave basis set 

A good example of orthonormal basis set, and one commonly employed in 
physics, is the plane-wave basis set. This basis set is closely related to Fourier 
transforms and it can be easily understood if concepts from Fourier analysis are 
known. 

A function f{x) defined on the entire real axis can be always expanded into 
Fourier components, f(k): 

f(x) = -= / f(k)e ikx dk (4.40) 



1 



oo 



f(k) = — / f(x)e-^dx. (4.41) 



For a function defined on a finite interval [—a/2, a/2], we can instead write 

f(x) = ~£/(fcn)^ (4-42) 



'a 

1 /-a/2 



f(k n ) = — / f(x)e" lk " x dx (4.43) 

"' -a/2 



where k n = 2im/a, n = 0, ±1, ±2, .... Note that the f(x) function of Eq.4.42 



is 



by construction a periodic function, with period equal to a: f(x + a) = f(x), as 
can be verified immediately. This implies that /(—a/2) = /(+a/2) must hold 
(also known under the name of periodic boundary conditions). The expressions 
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reported here are immediately generalized to three or more dimensions. In the 
following only a simple one-dimensional case will be shown. 



Let us define our plane-wave basis set bi(x) according to Eq.(4.42): 



1 2ir 

bi(x) = -=e lk * x , h = —i, i = 0,±l,±2,...,±N (4.44) 

a 

and the corresponding coefficients c$ for the wave function ip(x) as 

ra/2 

d= b*{x)ip{x)dx = (bi\ip), ip (x) = Y\cik(x). (4.45) 

J-a/2 ^ 

This base, composed of 2N + 1 functions, becomes a complete basis set in the 
limit N — > oo. This is a consequence of well-known properties of Fourier series. 
It is also straightforward to verify that the basis is orthonormal: Sij = {bt\bj) = 
5{j. The solution of the problem of a particle under a potential requires thus 
the diagonalization of the Hamiltonian matrix, whose matrix elements: 

2 

Hij = (bi\H\b 3 ) = + V(x)\bj) (4.46) 

can be trivially calculated. The kinetic term is diagonal (i.e. it can be repre- 
sented by a diagonal matrix): 

(bi\ — \bi) = / b*(x)-r4(x)dx = 5a -. (4.47) 

x '2m' 31 2m J - a /2 dx 2 y ' 3 2m v ; 

The potential term is nothing but the Fourier transform of the potential (apart 
from a multiplicative factor): 

(64|^(x)|6 7 -> = - r 2 V(x)e~ i{k ^ )x dx = -^=V{h - kA. (4.48) 

A known property of Fourier transform ensures that the matrix elements of the 
potential tend to zero for large values of ki — kj. The decay rate will depend 
upon the spatial variation of the potential: faster for slowly varying potentials, 
and vice versa. Potentials and wave functions varying on a typical length scale 
A have a significant Fourier transform up to k max ~ 2n/X. In this way we can 
estimate the number of plane waves needed to solve a problem. 



4.4 Code: pwell 

Let us consider the simple problem of a potential well with finite depth V$: 

V(x) = per z<-^,x>^ (4.49) 
V(x) = -V Q per - ^ < x < - (4.50) 
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with Vq > 0, b < a. The matrix elements of the Hamiltonian are given by 



Eq.(4.47) for the kinetic part, by Eq.(4.48) for the potential. The latter can be 



explicitly calculated: 



1 r b / 2 

(&i|F(x)|& 7 -) = — / Voe-^-^dx (4.51) 
a J-b/2 



a —i(ki — kj) 
V sin (b(ki - kj)/2) 



6/2 

V (4.52) 

-6/2 



a (ki — kj)/2 
The case ki = kj must be separately treated, yielding 



ki / kj. (4.53) 



V(0) = (4.54) 



a 



Code pwell.f9(j^] (or pwell.cj^]) generates the ki, fills the matrix Hij and di- 
agonalizes it. The code uses units in which h 2 /2m = 1 (e.g. atomic Rydberg 
units). Input data are: width (b) and depth (Vo) of the potential well, width 
of the box (a), number of plane waves (2N +1). On output, the code prints 
the three lowest energy levels; moreover it writes to file gs-wfc.out the wave 
function of the ground state. 

4.4.1 Diagonalization routines 



The practical solution of the secular equation, Eq.(4.32), is not done by naively 
calculating the determinant and finding its roots! Various well-established, 
robust and fast diagonalization algorithms are known. Typically they are based 
on the reduction of the original matrix to Hessenberg or tridiagonal form via 
successive transformations. All such algorithms require the entire matrix (or at 
least half, exploiting hermiticity) to be available in memory at the same time, 
plus some work arrays. The time spent in the diagonalization is practically 
independent on the content of the matrix and it is invariably of the order of 
0(N 3 ) floating-point operations for a N x N matrix, even if eigenvalues only 
and not eigenvectors are desired. Matrix diagonalization used to be a major 
bottleneck in computation, due to its memory and time requirements. With 
modern computers, diagonalization of 1000 x 1000 matrix is done in less than 
no time. Still, memory growing as ./V 2 and time as iV 3 are still a serious obstacle 
towards larger N. At the end of these lecture notes alternative approaches will 
be mentioned. 

The computer implementation of diagonalization algorithms is also rather 
well-established. In our code we use subroutine dsyev . ^] from the linear alge- 
bra library LAPACK^J Several subroutines from the basic linear algebra library 



1 http://www.fisica.uniud.it/%7Egiannozz/Corsi/MQ/Software/F90/pwell.f90 



http:/ /www.fisica.uniud.it/%7Egiarmozz/Corsi/MQ/Software/C/pwell.c 



http:/ /www.fisica.uniud.it/%7Egiannozz/Corsi/MQ/Software/Lapack/dsyev.f 



http://www.netlib.org/lapack/! 
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(collected here: dgemm . are also required, dsyev implements re- 
duction to tri-diagonal form for a real symmetric matrix (d=double precision, 
sy=symmetric matrix, ev=calculate eigenvalues and eigenvectors). The usage 
of dsyev requires either linking to a pre-compiled LAPACK and BLAS libraries, 
or compilation of the Fortran version and subsequent linking. Instructions on 
the correct way to call dsyev are contained in the header of the subroutine. 
Beware: most diagonalization routines overwrite the matrix! 

For the C version of the code, it may be necessary to add an underscore (as in 



dsyev_() ) in the calling program. Moreover, the tricks explained in Sec 0.1.4 
are used to define matrices and to pass arguments to BLAS and LAPACK 
routines. 

4.4.2 Laboratory 

• Observe how the results converge with respect to the number of plane 
waves, verify the form of the wave function. Verify the energy you get for 
a known case. You may use for instance the following case: for Vq = 1, 
6 = 2, the exact result is E = —0.4538. You may (and should) also verify 
the limit Vq — > oo (what are the energy levels?). 

• Observe how the results converge with respect to a. Note that for values 
of a not so big with respect to b, the energy calculated with the variational 
method is lower than the exact value. Why is it so? 

• Plot the ground-state wavefunction. You can also modify the code to 
write excited states. Do they look like what you expect? 

Possible code modifications and extensions: 

• Modify the code, adapting it to a potential well having a Gaussian form 
(whose Fourier transform can be analytically calculated: what is the 
Fourier transform of a Gaussian function?) For the same "width", which 
problem converges more quickly: the square well or the Gaussian well? 

• We know that for a symmetric potential, i.e. V(—x) = V(x), the solutions 
have a well-defined parity, alternating even and odd parity (ground state 
even, first excited state odd, and so on). Exploit this property to reduce 
the problem into two subproblems, one for even states and one for odd 
states. Use sine and cosine functions, obtained by suitable combinations 
of plane waves as above. Beware the correct normalization and the k n = 
term! Why is this convenient? What is gained in computational terms? 



"http:/ /www. netlib.org/blas/ 



http:/ /www.nsica.uniud.it/%7Egiannozz/Corsi/MQ/Software/Blas/dgemm.f 
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Chapter 5 

Non-orthonormal basis sets 



We consider here the extension of the variational method to the case of non- 
orthonormal basis sets. We consider in particular the cause of Gaussian func- 
tions as basis set. This kind of basis set is especially used and useful in molecular 
calculations using Quantum Chemistry approaches. 



5.1 Non-orthonormal basis set 

Linear-algebra methods allow to treat with no special difficulty also the case in 
which the basis is formed by functions that are not orthonormal, i.e. for which 



Sij = (bilbj) = J b*bjdv (5.1) 

is not simply equal to Sij. The quantities Sij are known as overlap integrals. 
It is sometimes practical to work with basis of this kind, rather than with an 
orthonormal basis. 

In principle, one could always obtain an orthonormal basis set from a non- 
orthonormal one using the Gram-Schmid orthogonalization procedure. An or- 
thogonal basis set bi is obtained as follows: 

h = h (5.2) 
b 2 = b 2 -bi(bi\b 2 )/(bi\bi) (5.3) 
h = 63-6 2 (6 2 |63>/(62|6 2 )-6 1 (6 1 |5 3 >/(6i|6i) (5.4) 

and so on. In practice, this procedure is seldom used and it is more convenient 
to follow a similar approach to that of Sect 4.2 In this case, Eq.(4.27) is 
generalized as 

G(d, ...,C N ) = c *i c j( H ij - eS ij) ( 5 - 5 ) 

and the minimum condition ( |4.31[ ) becomes 

Y,(Hij - eSij)cj = (5.6) 
j 

or, in matrix form, 

He = eSc (5.7) 
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known as generalized eigenvalue problem. 

The solution of a generalized eigenvalue problem is in practice equivalent to 
the solution of two simple eigenvalue problems. Let us first solve the auxiliary 
problem: 

Sd = ad (5.8) 



completely analogous to the problem (4.33). We can thus find a unitary matrix 
D (obtained by putting eigenvectors as columns side by side), such that D~ l SD 
is diagonal (D~ l = D*), and whose non-zero elements are the eigenvalues a. 
We find an equation similar to Eq. (14.39 ) : 



D*k SijDj n — 0~ n 5kn- 



(5.9) 



Note that all a n > 0: an overlap matrix is positive definite. In fact, 

(b n \b n ), \K) = ^2 D jn\bj) 



(5.10) 



and | b) is the rotated basis set in which S is diagonal. Note that a zero eigenvalue 
a means that the corresponding \b) has zero norm, i.e. one of the b functions is 
a linear combination of the other functions. In that case, the matrix is called 
singular and some matrix operations (e.g. inversion) are not well defined. 
Let us define now a second transformation matrix 



A 



D 



ij 



/a. 



We can write 



^ A* ik SjjA 



J I' 



6, 



kn 



(5.11) 



(5.12) 



(note that A is not unitary) or, in matrix form, A^SA = I. Let us now define 

c = Av (5.13) 



With this definition, Eq.(5.7) becomes 

HAv = eSAv 

We multiply to the left by A^: 

A^HAv = eA^SAv = ev 



(5.14) 



(5.15) 



Thus, by solving the secular problem for operator A^HA, we find the desired 
eigenvalues for the energy. In order to obtain the eigenvectors in the starting 
base, it is sufficient, following Eq.(5.13 ), to apply operator A to each eigenvector. 
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5.1.1 Gaussian functions 



Gaussian functions are frequently used as basis functions, especially for atomic 
and molecular calculations. They are known as GTO: Gaussian- Type Orbitals. 
An important feature of Gaussians is that the product of two Gaussian func- 
tions, centered at different centers, can be written as a single Gaussian: 

axi + /3r 2 



e _ a ( r __ ri )2 e _ /3(r _ r2 )2 



-(a+P)(r-r y 



(ri-r 2 y 



Some useful integrals involving Gaussian functions: 



f 

Jo 



dx 





/ 




Jo 



xe 



dx 



2a 



J o 



a + j3 



1 

2a' 



from which one derives 



e- ax2 x 2n dx 



e~ ax x 2n+1 dx 



da n jo 

f)n r oo 

j 



dx 



xe 



dx 



(2n-l)!!^/2 

2n+l Q ,n+l/2 

n! 



o 



2a n+1 



(5.16) 

(5.17) 

(5.18) 
(5.19) 



5.1.2 Exponentials 



Basis functions composed of Hydrogen-like wave functions (i.e. exponentials) 
are also used in Quantum Chemistry as alternatives to Gaussian functions. 
They are known as STO: Slater- Type Orbitals. Some useful orbitals involving 
STOs: 



-2Zr 

d 3 r = 4tt I r< 

r Jo 



~ 2Zr dr = 4tt 



-2Zr 



1 



-2Z(n+r 2 ) 



|ri - r 2 | 

5.2 Code: hydrogen_gauss 



-d i rid i r 2 



2Z 4Z 2 



5vr 2 



7T 



8Z 5 ' 



(5.20) 



(5.21) 



Code hydrogen.gauss . f9C0 (or hydrogen.gauss . cj^]) solves the secular prob- 
lem for the hydrogen atom using two different non-orthonormal basis sets: 



1. a Gaussian, "S-wave" basis set: 

bi(r) 



(5.22) 



2. a Gaussian "P-wave" basis set, existing in three different choices, corre- 
sponding to the different values m of the projection of the angular mo- 
mentum L z : 

b i {r) = xe-^ r \ b i {v)=ye-^ r \ 6,(r) = ze~^\ (5.23) 
(actually only the third choice corresponds to a well-defined value m = 0) 



1 http:/ /www. fisica.uniud.it/%7Egiannozz/Corsi/MQ/Software/F90/hydrogen_gauss.f90 



http:/ /www.fisica.uniud.it/%7Egiannozz/Corsi/MQ/Software/C/hydrogen_gauss.c 
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The Hamiltonian operator for this problem is obviously 



H 



h 2 V 2 Zq 2 e 
2m P r 



(5.24) 



For the hydrogen atom, Z = 1. 

Calculations for S- and P-wave gaussians are completely independent. In 
fact, the two sets of basis functions are mutually orthogonal: Sij = if i is a 
S-wave, j is a P-wave gaussian, as evident from the different parity of the two 
sets of functions. Moreover the matrix elements Hij of the Hamiltonian are 
also zero between states of different angular momentum, for obvious symmetry 
reasons. The S and H matrices are thus block matrices and the eigenvalue 
problem can be solved separately for each block. The P-wave basis is clearly 
unfit to describe the ground state, since it doesn't have the correct symmetry, 
and it is included mainly as an example. 

The code reads from file a list of exponents, ai, and proceeds to evaluate all 
matrix elements and Sij. The calculation is based upon analytical results 



for integrals of Gaussian functions (Sec 5.1.1 ). In particular, for S-wave one has 

3/2 



s, 



IT 



ai + aj 



(5.25) 



while the kinetic and Coulomb terms in Hij are respectively 



H 



K 





h 2 v 2 ' 


e 




h 2 Qaiaj 


2m e 


2m e ai + aj 




-1- 




\ Zqf 


e -^ d 3 r = _ 




r 



IT 



3/2 



ai + a,- 



2-irZq] 



ai + a 



(5.26) 



(5.27) 
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For the P-wave basis the procedure is analogous, using the corresponding (and 
more complex) analytical expressions for integrals. 

The code then calls subroutine diag that solves the generalized secular 
problem (i.e. it applies the variational principle). Subroutine diag returns a 
vector e containing eigenvalues (in order of increasing energy) and a matrix v 
containing the eigenvectors, i.e. the expansion coefficients of wave functions. 

Internally, diag performs the calculation described in the preceding section 
in two stages. The solution of the simple eigenvalue problem is performed by 



the subroutine dsyev we have already seen in Sect. 4.4 

In principle, one could use a single LAPACK routine, dsygv, that solves 
the generalized secular problem, Hip = eSip, with a single call. In practice, 
one has to be careful to avoid numerical instabilities related to the problem 



of linear dependencies among basis functions (see Eq.(5.10) and the following 
discussion). Inside routine diag, all eigenvectors of matrix S corresponding to 
very small eigenvectors, i.e. smaller than a pre-fixed threshold are thrown away, 
before proceeding with the second diagonalization. The number of linearly 
independent eigenvectors is reprinted in the output. 

The reason for such procedure is that it is not uncommon to discover that 
some basis functions can almost exactly be written as sums of some other basis 
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functions. This does not happen if the basis set is well chosen, but it can happen 
if the basis set functions are too many or not well chosen (e.g. with too close 
exponents). A wise choice of the ay coefficients is needed in order to have a 
high accuracy without numerical instabilities. 

The code then proceeds and writes to files s-coef f . out (or p-coef f . out) 
the coefficients of the expansion of the wave function into Gaussians. The 
ground state wave function is written into file s-wfc.out (or p-wfc.out). 

Notice the usage of dgemm calls to perform matrix-matrix multiplication. 
The header of dgemm. f contains a detailed documentation on how to call it. 
Its usage may seem awkward at a first (and also at a second) sight. This 
is a consequence in part of the Fortran way to store matrices (requiring the 
knowledge of the first, or "leading", dimension of matrices); in part, of the 
old-style Fortran way to pass variables to subroutines only under the form of 
pointers (also for vectors and arrays). Since the Fortran-90 syntax and MATMUL 
intrinsic are so much easier to use: why bother with dgemm and its obscure 
syntax? The reason is efficiency: very efficient implementations of dgemm exist 
for modern architectures. 

For the C version of the code, and how matrices are introduced and passed 



to Fortran routines, see Sec 0.1.4 



5.2.1 Laboratory 

• Verify the accuracy of the energy eigenvalues, starting with 1 Gaussian, 
then 2, then 3. Try to find the best values for the coefficients for the Is 
state (i.e. the values that yield the lowest energy). 

• Compare with the the solutions obtained using code hydrogen_radial. 
Plot the Is numerical solution (calculated with high accuracy) and the 
"best" Is solution for 1, 2, 3, Gaussians (you will need to multiply the 
latter by a factor \f^n: why? where does it come from?). What do you 
observe? where is the most significant error concentrated? 

• Compare with the results for the following optimized basis sets (a.u.): 

- three Gaussians: a\ = 0.109818, a 2 = 0.405771, a 3 = 2.22776 
(known as "STO-3G" in Quantum-Chemistry jargon) 

- four Gaussians: a x = 0.121949, a 2 = 0.444529, a 3 = 1.962079, 
a 4 = 13.00773 

• Observe and discuss the ground state obtained using the P-wave basis set 

• Observe the effects related to the number of basis functions, and to the 
choice of the parameters a. Try for instance to choose the characteristic 
Gaussian widths, A = l/^/a, as uniformly distributed between suitably 
chosen X min and X ma x- 

• For Z > 1 , how would you re-scale the coefficients of the optimized Gaus- 
sians above? 
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Chapter 6 

Self-consistent Field 



A way to solve a system of many electrons is to consider each electron under the 
electrostatic field generated by all other electrons. The many-body problem is 
thus reduced to the solution of single-electron Schrodinger equations under an 
effective potential. The latter is generated by the charge distribution of all other 
electrons in a self- consistent way. This idea is formalized in a rigorous way in 
the Hartree-Fock method and in Density-Functional theory. In the following 
we will see an historical approach of this kind: the Hartree method. 



6.1 The Hartree Approximation 

The idea of the Hartree method is to try to approximate the wave functions, 
solution of the Schrodinger equation for N electrons, as a product of single- 
electron wave functions, called atomic orbitals. As we have seen, the best 
possible approximation consists in applying the variational principle, by mini- 
mizing the expectation value of the energy E = (ip\H\ip) for state \ip). 

The Hamiltonian of an atom having a nucleus with charge Z and N electrons 

is 

where the sum is over pairs of electrons i and j, i.e. each pair appears only 
once. Alternatively: 

N-l N 

E = E E (6-2) 

{ij) i=l j=i+l 

It is convenient to introduce one-electron and two-electrons operators: 

fi - ( 6 - 3 ) 
2m e ri 

9i3 - f (6-4) 
'ij 

With such notation, the Hamiltonian is written as 

# = E/* + E^ (6-5) 

* (ij) 
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provided that g act symmetrically on the two electrons (definitely true for the 
Coulomb interaction). 



6.2 Hartree Equations 

Let us assume that the total wave function can be expressed as a product of 
single-electron orbitals (assumed to be orthonormal) : 



^(l,2,...,iV-) = 1 (l)0 2 (2)...0 A r(iV) 



(6.6) 
(6.7) 



Variables 1, 2, .., mean position and spin variables for electrons 1,2,...; J dvi 
means integration on coordinates and sum over spin components. Index j labels 
instead the quantum numbers used to classify a given single-electron orbitals. 
All orbitals must be different: Pauli's exclusion principle tells us that we cannot 
build a wave function for many electrons using twice the same single-electron 
orbital. In practice, orbitals for the case of atoms are classified using the main 
quantum number n, orbital angular momentum t and its projection m. 
The expectation value of the energy is 



\H\ 



. . .<p N (N) dv 1 ...dv N 



dvidvj 



j #(l)/i&(l)dt>i + X; / 0*(l)^(2) 9l2 ^(l)^(2)^ 1 ^ 2 

* (*i> 

(6.8) 



In the first step, we made use of orthonormality (6.7); in the second we just 



renamed dummy integration variables with the "standard" choice 1 and 2. 



Let us now apply the variational principle to formulation (4.10), with the 
constraints that all integrals 



h = j 4>m<t>k(l)dvi (6.9) 
are constant, i.e. the normalization of each orbital function is preserved: 



- 5> fc I fc ) =0 
k J 



(6.10) 



where e& are Lagrange multipliers, to be determined. Let us vary only the 
orbital function <bu. We find 



5I k = / 6<f>%(l)<f> k (i) dvi + c.c. 



(6.11) 
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(the variations of all other normalization integers will be zero) and, using the 
hermiticity of H as in Sec 4.1.1 



8(i>\H\i>) 



<W)/i<^(l)cfoi + c.c. 



+ I ^(1)^(2)512^(1)0,(2) ^i^ 2 + c.c. (6.12) 



This result is obtained by noticing that the only involved terms of Eq.(6.8) are 



those with i = k or j = k, and that each pair is counted only once. For instance, 
for 4 electrons the pairs are 12, 13, 14, 23, 24, 34; if I choose k = 3 the only 
contributions come from 13, 23, 34, i.e. J2j^=k (since g is a symmetric operator, 
the order of indices in a pair is irrelevant) 

Thus the variational principle takes the form 



j^k 



dv\ + c.c. 







i.e. the term between square brackets (and its complex conjugate) must vanish 
so that the above equation is satisfied for any variation: 



h<(> k (l) + E / 4> -(2)512^(1)^(2) dv 2 = e k </> h (l) 



.13) 



These are the Hartree equations (k = 1,...,N). It is useful to write down 
explicitly the operators: 



1,2 vlMV - —Mi) + 



2m, 



E 



b*(2)^{2)dv 2 



(1) =e*0fc(l) 



(6.14) 

We remark that each of them looks like a Schrodinger equation in which in 
addition to the Coulomb potential there is a Hartree potential: 



Vff(ri) = / Pk (2)^-dv2 
r\2 



(6.15) 
(6.16) 



where we have used 

Pfc (2) = E^(2)^(2) 
j^k 

Pj is the charge density due to all electrons differing from the one for which we 
are writing the equation. 



6.2.1 Eigenvalues and Hartree energy 



Let us multiply Hartree equation, Eq(6.13), by <ft k (l), integrate and sum over 
orbit als: we obtain 



E £ fc = E / ^(i)/i0 fc (i)^i + EE / ^(1)^(1)512^(2)^(2)^^2. 

(6.17) 



k 3+k' 
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Let us compare this expression with the energy for the many-electron system, 



Eq.(6.8). The Coulomb repulsion energy is counted twice, since each < jk > 
pair is present twice in the sum. The energies thus given by the sum of eigen- 
values of the Hartree equation, minus the Coulomb repulsive energy: 

S = E e *- E f ^k(l)M^9i2<f>*(2)4> j ( 2 )dndv 2 . (6.18) 
k <jk> J 

6.3 Self-consistent potential 



Eq.(6.15) represents the electrostatic potential at point ri generated by a charge 
distribution This fact clarifies the meaning of the Hartree approximation. 
Assuming that ip is factorizable into a product, we have formally assumed that 
the electrons are independent. This is of course not true at all: the electrons 
are strongly interacting particles. The approximation is however not so bad if 
the Coulomb repulsion between electrons is accounted for under the form of an 
average field Vh , containing the combined repulsion from all other electrons on 
the electron that we are considering. Such effect adds to the Coulomb attraction 
exerted by the nucleus, and partially screens it. The electrons behave as i/they 
were independent, but under a potential —Zq^/r + Vn(r) instead of —Zq\jr of 
the nucleus alone. 

Vh{t) is however not a "true" potential, since its definition depends upon 
the charge density distributions of the electrons, that depend in turn upon the 
solutions of our equations. The potential is thus not known a priori, but it is a 
function of the solution. This type of equations is known as integro-differential 
equations. 

The equations can be solved in an iterative way, after an initial guess of the 
orbitals is assumed. The procedure is as follows: 

1. calculate the charge density (sum of the square moduli of the orbitals) 

2. calculate the Hartree potential generated by such charge density (using 
classical electrostatics) 

3. solve the equations to get new orbitals. 

The solution of the equations can be found using the methods presented in 
Chj2j The electron density is build by filling the orbitals in order of increasing 
energy (following Pauli's principle) until all electrons are "placed". 

In general, the final orbitals will differ from the starting ones. The procedure 
is then iterated - by using as new starting functions the final functions, or with 
more sophisticated methods - until the final and starting orbitals are the same 
(within some suitably defined numerical threshold). The resulting potential Vh 
is then consistent with the orbitals that generate it, and it is for this reason 
called self- consistent field. 

6.3.1 Self-consistent potential in atoms 

For closed-shell atoms, a big simplification can be achieved: Vh is a central 
field, i.e. it depends only on the distance r\ between the electron and the 
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nucleus. Even in open-shell atoms, this can be imposed as an approximation, 
by spherically averaging p^. The simplification is considerable, since we know 



a priori that the orbitals will be factorized as in Eq.(2.9). The angular part is 

and m, while 
Of course the 



given by spherical harmonics, labelled with quantum numbers 
the radial part is characterized by quantum numbers n and t 
accidental degeneracy for different 0. is no longer present. It turns out that even 
in open-shell atoms, this is an excellent approximation. 

Let us consider the case of two-electron atoms. The Hartree equations, 
Eq.(6.14), for orbital k = 1 reduces to 



ti 2 



2m, 



-v?<Mi) 



zfi 



01 (1) + 



ri2 



-M2)dv 2 



01 (1) = ei0i(l) (6.19) 



For the ground state of He, we can assume that (j>i and 4>2 have the same spher- 
ically symmetric coordinate part, (j)(r), and opposite spins: <p\ = 4>(r)v + (a), 
4>2 = </>(r)v—(cr). Eq.(6.19) further simplifies to: 



2m, 



V?0(n) 



9^1 

ri2 



(r 2 )| 2 d 3 r 2 



0(n) = e0(n) (6.20) 



6.4 Code: helium_hf_radial 

Code heliumJif _radial . f 9CQ (or heliumJif _radial . cj^]) solves Hartree equa- 
tions for the ground state of He atom. helium_hf .radial is based on code 
hydrogen_radial and uses the same integration algorithm based on Numerov's 
method. The new part is the implementation of the method of self-consistent 
field for finding the orbitals. 

The calculation consists in solving the radial part of the Hartree equation 



(6.20). The effective potential V sc f is the sum of the Coulomb potential of the 



nucleus, plus the (spherically symmetric) Hartree potential 

Vscf{r) = + V H (r), V„(n) = q 2 e f (6.21) 

r J r\2 

We start from an initial estimate for Vn(r), calculated in routine init.pot ( 
vjj\r) = 0, simply). With the ground state R(r) obtained from such potential, 
we calculate in routine rho_of r the charge density p(r) = |i?(r)| 2 /47r (note that 
p is here only the contribution of the other orbital, so half the total charge, and 
remember the presence of the angular part!). Routine v_of _rho re-calculates 
the new Hartree potential (r) by integration, using the Gauss theorem: 

VZ u \r) = V + q 2 e f ^fids, Q(s) = [ p(r)4irr 2 dr (6.22) 

Jr ma x s Jr<s 

where Q(s) is the charge contained in the sphere of radius s; r max is the out- 
ermost grid point, such that the potential has already assumed the asymptotic 
value Vq = qe/r max , valid at large r. 



1 http://www.fisica.uniud.it/%7Egiannozz/Corsi/MQ/Software/F90/helium_hf_radial.f90 



http:/ /www.fisica.uniud.it/%7Egiannozz/Corsi/MQ/Software/C/helium_hf_radial.c 
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The Hartree potential is then re-introduced in the calculation not directly 
but as a linear combination of the old and the new potential. This is the 
simplest technique to ensure that the self- consistent procedure converges. It 
is not needed in this case, but in most cases it is: there is no guarantee that 
re-inserting the new potential in input will lead to convergence. We can write 

yvn,ne W{r) = ^yout (r) + (1 _ ^(r), ( 6 .23) 

where < /3 < 1. The procedure is iterated (repeated) until convergence is 
achieved. The latter is verified on the norm, i.e. the integral of the square, of 
VjP*(r) — Vjp(r). square. 

In output the code prints the eigenvalue ei of Eq |6.13 plus various terms of 



the energy, with rather obvious meaning except the term Variational correction. 
This is 

5E= !{V£{r) - V£ ut (r))p(r)d?r (6.24) 



and it is useful to correctj^] the value of the energy obtained by summing the 
eigenvalues as in Eq.( 6.18| , so that it is the same as the one obtained using 



Eq.(6.8), where eigenvalues are not used. The two values of the energy are 
printed side by side. Also noticeable is the "virial check": for a Coulomb 
potential, the virial theorem states that (T) = — (V)/2, where the two terms 
are respectively the average values of the kinetic and the potential energy. It 
can be demonstrated that the Hartree-Fock solution obeys the virial theorem. 

6.4.1 Laboratory 

• Observe the behavior of self-consistency, verify that the energy (but not 
the single terms of it!) decreases monotonically. 

• Compare the energy obtained with this and with other methods: per- 
turbation theory with hydrogen-like wave functions (E = —5.5 Ry, Sect. 



D.l), variational theory with effective Z (E = —5.695 Ry, Sect. D.3), 
best numerical Hartree (-Fock) result (E = —5.72336 Ry, as found in the 
literature), experimental result (E = —5.8074 Ry). 

Make a plot of orbitals (file wfc.out) for different n and i. Note that 
the orbitals and the corresponding eigenvalues become more and more 
hydrogen-like for higher n. Can you explain why? 

If you do not know how to answer the previous question: make a plot 
of V sc j (file pot. out) and compare its behavior with that of the —Zq\jr 
term. What do you notice? 

Plot the Is orbital together with those calculated by hydrogen_radial 



for Hydrogen [Z = 1), He+ (Z = 2), and for a Z = 1.6875. See Sect. |D.3 
if you cannot make sense of the results. 



Eigenvalues are calculated using the input potential; the other terms are calculated using 
output potential 
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Chapter 7 

The Hartree-Fock 
approximation 



The Hartree method is useful as an introduction to the solution of many-particle 
system and to the concepts of self-consistency and of the self-consistent-field, 
but its importance is confined to the history of physics. In fact the Hartree 
method is not just approximate: it is wrong, by construction, since its wave 
function is not antisymmetric! A better approach, that correctly takes into 
account the antisymmetric character of the the wave functions is the Hartree- 
Fock approach. The price to pay is the presence in the equations of a non local, 
and thus more complex, exchange potential. 



7.1 Hartree-Fock method 



Let us re-consider the Hartree wave function. The simple product: 
M1,2,...,N) = &(l)fo(2)... <j> N (N) 



(7.1) 



does not satisfy the principle of indistinguishability, because it is not an eigen- 
state of permutation operators. It is however possible to build an antisymmetric 
solution by introducing the following Slater determinant: 



V>(l,...,Af) 



1 



Mi) 



<Mi) 



(7.2) 



The exchange of two particles is equivalent to the exchange of two columns, 
which induces, due to the known properties of determinants, a change of sign. 
Note that if two rows are equal, the determinant is zero: all (f>iS must be 
different. This demonstrates Pauli's exclusion principle, two (or more) identical 
fermions cannot occupy the same state. 

Note that the single-electron orbitals 4>i are assumed to be orthonormal: 



J #(1)^(1)^! = ^ 



(7.3) 
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where the "integral" on dv\ means as usual "integration on coordinates, sum 
over spin components" . We follow the same path of Sec. (6.2) used to derive 
Hartree equations, Eq.(6.13 ). Since a determinant for N electrons has AH terms, 
we need a way to write matrix elements between determinants on a finite paper 
surface. The following property, valid for any (symmetric) operator F and 
determinantal functions ip and ij/, is very useful: 



i 

m 



0i (i) • Fi(N) 



<j>* N (l) . <p* N (N) 

<Pl(l)...<P* N (N)F' 



#v(l) ■ 
#(JV) 



^(1) . <^(iV) 



01 (iV) 
^v(^) 
dui . . . dvN 



dvi . . . dvj\[ 



(7.4) 



(by expanding the first determinant, one gets N\ terms that, once integrated, 
are identical). From the above property it is immediate (and boring) to obtain 
the matrix elements for one- and two-electron operators: 



(VI = £ / #(i)/i&(i)<foi 

i i 

(as in the Hartree approximation), and 



(7.5) 



(ij) 



E 

(ij) 



#(l)#(2) 5 i 2 [&(1)&(2) - ^-(1)^(2)] d«idu 2 (7.6) 



The integrals implicitly include summation over spin components. If we assume 
that g\2 depends only upon coordinates (as in Coulomb interaction) and not 
upon spin, the second term: 



<j>* (!)<£* (2)012^- (1)^(2) dmd«2 



(7.7) 



is zero is i and j states are different (the spin parts are not affected by g±2 and 
they are orthogonal if relative to different spins) . 

It is convenient to move to a scheme in which the spin variables are not 
explicitly included in the orbital index. Eq.(7.6) can then be written as 

MI>«M=E / 0*(1)0K 2 )S12 [&Wi(2) -5(^,^)^(1)^(2)] d Vl dv2 

(ij) (ij) 

(7.8) 

where Uj is the spin of electron i, and: 

5(ai,aj) = if Ui 7^ crj 



1 if ov 



(Jo- 



in summary: 



]T | 0?(l)/ a 0,-(l)d«i 



(7.9) 



+ E 



^(l)^(2) 5 i 2 [^(1)^(2) - 5(^,^)^(1)^(2)] dv!dv2 
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Now that we have the expectation value of the energy, we can apply the vari- 
ational principle. In principle we must impose normalization constraints such 
that not only all <$>% stay normalized (as we did in the derivation of Hartree's 
equation) but also all pairs (pi, with same spin are orthogonal, i.e., a (tri- 
angular) matrix e^- of Lagrange multipliers would be needed. It can be shown 
however (details e.g. on Slater's book, Quantum theory of matter) that it is 
always possible to find a transformation to a solution in which the matrix of 
Lagrange multipliers is diagonal. We assume that we re dealing with such a 
case. 

Let us omit the details of the derivation and move to the final Hartree-Fock 
equations: 



MkiX) + E / <^( 2 )3i2 fe(l)^-(2)-%,a J )^(l)^(2)] dv 2 = e k <t> k (l) 

(7.10) 



3 

or, in explicit form, 

h 2 2 Zq* ^ f qf 



-V? fail) - + (7-11) 

- 5(a k ,o-j)(/>k(2)(t)j(l)] dv 2 = e k (f> k (l) 



The energy of the system, Eq. 7.9, can be expressed, analogously to the Hartree 



case, via the sum of eigenvectors of Eq.(7.11), minus a term compensating the 



double counting of Coulomb repulsion and of the exchange energy: 
£ = E £ fc- E / #fe(l)$(2)si2 IMVM 2 ) - 8{a h a k )cj )j {l) ( t> k (2)}dv x dv 2 . 

k <jk> 

(7-12) 

Eq.(7.11 ) has normally an infinite number of solutions, of which only the lowest- 
energy N will be occupied by electrons, the rest playing the role of excited 
states. The sum over index j runs only on occupied states. 

Let us carefully observe the differences with respect to Hartree equations, 



Eqs.(6.13): 



1. J2j a l so includes the case j = k. 

2. for electrons j with same spin as that of k, there is an additional exchange 
term 

3. due to the exchange term, the case j = k yields a nonzero contribution 
only if the spins are different. 

7.1.1 Coulomb and exchange potentials 

Let us analyze the physical meaning of the Hartree-Fock equations. We re- write 
them under the form 

1,2 -Vfr fc (l) - ^W) + VkWfc(l) + (VMO) = (7-13) 



2m e n 
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where we have introduced a "Hartree potential" Vh (it is not the same as in the 
Hartree equations!) and an "exchange potential" V x . The Hartree potential is 
the same for all orbitals: 

2 2 

Vff(l) = E / ^(2)^^(2)^ 2 = / p(2)^dv 2 , (7.14) 

where we have introduce the charge density 

p(2) = ]>>*(2)^(2). (7.15) 

j 

We can verify that p(l) is equal to the probability to find an electron in (1), 
that is, 

p(l) = N [ |*(1, 2, .., iV)! 2 ^...^. (7.16) 



The exchange term: 

(V x <f> k )(l) = -5>fa.o-i) / ^(2)^^(2)^(1) d«a (7.17) 

does not have the simple form of the Hartree potential: Vfl-(l)</>fc(l), where 
Vff(l) comes from an integration over variable 2. It has instead a form like 



{V x <j>){l) = I ^(1,2)^(2)^2 (7.18) 
typical of a non local interaction. 



7.1.2 Correlation energy 

The Hartree- Fock solution is not exact: it would be if the system under study 
were described by a wave function having the form of a Slater determinant. This 
is in general not true. The energy difference between the exact and Hartree- 
Fock solution is known as correlation energy^ The origin of the name comes 
from the fact that the Hartree- Fock approximation misses part of the "electron 
correlation" : the effects of an electron on all others. This is present in Hartree- 
Fock via the exchange and electrostatic interactions; more subtle effects are not 
accounted for, because they require a more general form of the wave function. 
For instance, the probability P(ri,r 2 ) to find an electron at distance r\ and 
one as distance r 2 from the origin is not simply equal to p(r\)p{r2), because 
electrons try to "avoid" each other. The correlation energy in the case of He 
atom is about 0.084 Ry: a small quantity relative to the energy (~ 1.5%), but 
not negligible. 

An obvious way to improve upon the Hartree- Fock results consists in allow- 
ing contributions from other Slater determinants to the wave function. This is 
the essence of the "configuration interaction" (CI) method. Its practical ap- 
plication requires a sophisticated "technology" to choose among the enormous 

1 Feynman called it stupidity energy, because the only physical quantity that it measures is 
our inability to find the exact solution! 
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number of possible Slater determinants a subset of most significant ones. Such 
technique, computationally very heavy, is used in quantum chemistry to get 
high-precision results in small molecules. Other, less heavy methods (the so- 
called M0ller-Plesset, MP, approaches) rely on perturbation theory to yield a 
rather good estimate of correlation energy. A completely different approach, 
which produces equations that are reminiscent of Hartree-Fock equations, is 
Density- Functional Theory (DFT), much used in condensed-matter physics. 



7.1.3 The Helium atom 



The solution of Hartree-Fock equations in atoms also commonly uses the cen- 



tral field approximation. This allows to factorize Eqs.(7.11) into a radial and 



an angular part, and to classify the solution with the "traditional" quantum 
numbers n, t, m. 



For He, the Hartree-Fock equations, Eq.(7.11), reduce to 
-.2 



2m, 



-V?0 X (1) - ^0i(l) + / [0i(2)0i(l) - 0i(2)0i(l)] dv 2 

n J 7*12 



+ / 05(2)— [02(2)0i (1) - S(ax, (72)0! (2)02(1)] dvz = d0i(l) (7.19) 

J r 12 

Since the integrand in the first integral is zero, 



2m, 



n 



2 

0i(l)+ / 0£(2)-^[0 2 (2)0 1 (l) 



-<5(<7i, <t 2 )0i (2)02(1)] du 2 = ei0i(l). 



(7.20) 



In the ground state, the two electrons have opposite spin (5(a±,a2) = 0) and 
occupy the same spherically symmetric orbital (that is: 0i and 02 are the same 
function, 0). This means that the Hartree-Fock equation, Eq.(7.20), for the 



ground state is the same as the Hartree equation, Eq.(6.19). In fact, the two 



electrons have opposite spin and there is thus no exchange. 

In general, one speaks of Restricted Hartree-Fock (RHF) for the frequent 
case in which all orbitals are present in pairs, formed by a same function of r, 
multiplied by spin functions of opposite spin. In the following, this will always 
be the case. 



7.2 Code: helium_hf_gauss 

The radial solution of Hartree-Fock equations is possible only for atoms or in 
some model systems. In most cases the solution is found by expanding on a 
suitable basis set, in analogy with the variational method. 

Let us re-write the Hartree-Fock equation - for the restricted case, i.e. no 
total spin - in the following form: 

^0 fc = e fc fc , k = l,...,N/2 (7.21) 

The index k labels the coordinate parts of the orbitals; for each k there is 
a spin-up and a spin-down orbital. T is called the Fock operator. It is of 
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course a non-local operator which depends upon the orbitals <pk- Let us look 
now for a solution under the form of an expansion on a basis of functions: 
</>fc(r) = J2i c | "«( r )- We find the Rothaan-Hartree-Fock equations: 

Fc^ = e k Sc^ (7.22) 

where = (c^ , ■ ■ ■ , cff) is the vector of the expansion coefficients, S is 
the superposition matrix, F is the matrix of the Fock operator on the basis set 
functions: 



Fij = {bi\F\bj), 
that after some algebra can be written as 

N/2 



S; 



*fc = /*+££ US 



S k >Jk) 



9iljr 



. 9ijlm 



I m \ k=l 

where, with the notations introduced in this chapter: 

j3 , 



1 1. 1 



9iljr 



b*(r 1 )b j {r 1 )g 1 2^(r2)b m (r 2 )d 3 r 1 d i r2. 



(7.23) 
(7.24) 

(7.25) 
(7.26) 



The sum over states between parentheses in Eq.(7.24) is called density matrix. 



The two terms in the second parentheses come respectively from the Hartree 
and the exchange potentials. 



The problem of Eq.(7.22) is more complex than a normal secular problem 



solvable by diagonalization, since the Fock matrix, Eq.(7.24), depends upon its 



own eigenvectors. It is however possible to reconduct the solution to a self- 
consistent procedure, in which at each step a fixed matrix is diagonalized (or, 
for a non-orthonormal basis, a generalized diagonalization is performed at each 
step). 

Code heliumJif _gauss . (or heliumJif _gauss . cj^]) solves Hartree-Fock 
equations for the ground state of He atom, using a basis set of S Gaussians. The 
basic ingredients are the same as in code hydrogen_gauss (for the calculation 
of single-electron matrix elements and for matrix diagonalization). Moreover 



we need an expression for the guj m matrix elements introduced in Eq.(7.26). 
Using the properties of products of Gaussians, Eq.(5.16), these can be written 
in terms of the integral 



-ii r: „— j3r& \ 



ri2 



d 6 nd 6 r2. 



(7.27) 



Let us look for a variable change that makes (ri — Y2) 2 to appear in the exponent 
of the gaussians: 



7 (ri - r 2 ) 2 + (ar 1 + 6r 2 ) 2 

a/3 



a + /3 



(ri -r 2 y + 




(7.28) 
(7.29) 



http:/ /www. fisica.uniud.it/%7Egiannozz/Corsi/MQ/Software/F90/helium_hLgauss.f90 



http:/ /www.fisica.uniud.it/%7Egiannozz/Corsi/MQ/Software/C/helium_hf_gauss.c 
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Let us now introduce a further variable change from (ri,r2) to (r, s), where 



The integral becomes 
/ = 



ri - r 2 , 



__a§_ r 2 1 a/3 „2 

g OL+0 _g CK + /3 

r 




9(ri,r 2 ) 



d 3 rd 3 s, 



(7.30) 



(7.31) 



<9(r,s) 

where the Jacobian is easily calculated as the determinant of the transformation 

3 



matrix, Eq.(7.30): 



<9(ri,r 2 ) 



d(r,s) 



\/a(3 
a + ft 



(7.32) 



The calculation of the integral is trivial and provides the required result: 

9iljr 



a 



(7.33) 



where a = on + aj, (3 = on + a m . 

The self-consistent procedure is even simpler than in code heliumJif _radial: 
at each step, the Fock matrix is re-calculated using the density matrix at the 
preceding step, with no special trick or algorithm, until energy converges within 
a given numerical threshold. 

7.2.1 Laboratory 

• Observe how the ground-state energy chang function of the number 

of Gaussians and of their coefficients. You may take the energy given 
by heliumJif _radial as the reference. Start with the 3 or 4 Gaussian 
basis set used for Hydrogen, then the same set of Gaussians with rescaled 
coefficients, i.e. such that they fit a Slater Is orbital for He + (Z = 2) and 



for a Hydrogen-like atom with Z = 1.6875 (see Sect. D.3) 



• Write to file the Is orbital (borrow the relative code from hydrogen_gauss) , 
plot it. Compare it with 

— the Is Slater orbital for Z = 1, Z = 1.6875, Z = 2, and 

— the Is orbital from code heliumJif .radial. 

Beware: there is a factor between the orbitals calculated with radial in- 
tegration and with a Gaussian basis set: which factor and why? 

• Try the following optimized basis set with four Gaussians, with coeffi- 
cients: 

ax = 0.297104, a 2 = 1.236745, q 3 = 5.749982, a 4 = 38.216677 a.u.. 
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Chapter 8 

Molecules 



8.1 Born-Oppenheimer approximation 

Let us consider a system of interacting nuclei and electrons. In general, the 
Hamiltonian of the system will depend upon all nuclear coordinates, R„, and 
all electronic coordinates, T{. For a system of n electrons under the field of iV 
nuclei with charge Z^, in principle one has to solve the following Schrodinger 
equation: 



where T/ is the kinetic energy of nuclei, Vn is the Coulomb repulsion between 
nuclei, V e j is the Coulomb attraction between nuclei and electrons, T e is the 
kinetic energy of electrons, V ee is the Coulomb repulsion between electrons: 



This looks like an impressive problem. It is however possible to exploit the mass 
difference between electrons and nuclei to separate the global problem into an 
electronic problem for fixed nuclei and a nuclear problem under an effective 
potential generated by electrons. Such separation is known as adiabatic or Born- 
Oppenheimer approximation. The crucial point is that the electronic motion 
is much faster than the nuclear motion: while forces on nuclei and electrons 
have the same order of magnitude, an electron is at least ~ 2000 times lighter 
than any nucleus. We can thus assume that at any time the electrons "follow" 
the nuclear motion, while the nuclei at any time "feel" an affective potential 
generated by electrons. Formally, we assume a wave function of the form 



where the electronic wave function Vr, ( r i) solves the following Schrodinger 
equation: 



(Tj + V n + V eI + T e + V ee ) *(R M , r*) = £^(R M , r t ) 



(8.1) 




*(R M ,r i ) = *(R„X ; (r i ) 



(8.3) 



(T e + Vee + V eI ) 1$ (n) = E^^(n). 



(8.4) 
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The index R is a reminder that both the wave function and the energy depend 
upon the nuclear coordinates, via V e j; the index I classifies electronic states. 
We now insert the wave function, Eq.(8.3), into Eq.(8.2) and notice that T e 
does not act on nuclear variables. We will get the following equation: 

[Tj + Vn + E®) *(R M )^( ri ) = EViR^gin). (8.5) 



If we now neglect the dependency upon R of the electronic wave functions in 
the kinetic term: 

Tj ($(R M )tf( ri )) * ^(rt) (2>*(R„)) . (8.6) 
we obtain a Schrodinger equation for nuclear coordinates only: 

(T I + V II + E^)^(K ll ) = E^(K ll ), (8.7) 

where electrons have "disappeared" into the eigenvalue E^ . The term Vjj+E^ 
plays the role of effective interaction potential between nuclei. Of course such 
potential, as well as eigenfunctions and eigenvalues of the nuclear problem, 
depends upon the particular electronic state. 

The Born-Oppenheimer approximation is very well verified, except the spe- 
cial cases of non-adiabatic phenomena (that are very important, though). The 
main neglected term in Eq 8.6 has the form 

and may if needed be added perturbation. 



8.2 Potential Energy Surface 

The Born-Oppenheimer approximation allows us to separately solve a Schro- 
dinger equation for electrons, Eq.(8.4), as a function of atomic positions, and 
a problem for nuclei only, Eq.(8.7). The latter is in fact a Schrodinger equa- 
tion in which the nuclei interact via an effective interatomic potential, V(R^) = 
Vn+E^\ a function of the atomic positions R^ and of the electronic state. The 
interatomic potential V(R /X ) is also known as potential energy surface ("poten- 
tial" and "potential energy" are in this context synonyms), or PES. It is clear 
that the nuclear motion is completely determined by the PES (assuming that 
the electron states does not change with time) since forces acting on nuclei are 
nothing but the gradient of the PES: 

F M = -V M F(R M ), (8.9) 

while equilibrium positions for nuclei, labeled with R^°\ are characterized by 
zero gradient of the PES (and thus of any force on nuclei): 

Ffl = -V„y(R<°>) = 0. (8.10) 
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In general, there can be many equilibrium points, either stable (a minimum: 
any displacement from the equilibrium point produces forces opposed to the 
displacement, i.e. the second derivative is positive everywhere) or unstable (a 
maximum or a saddle point: for at least some directions of displacement from 
equilibrium, there are forces in the direction of the displacement, i.e. there 
are negative second derivatives). Among the various minima, there will be a 
global minimum, the lowest-energy one, corresponding to the ground state of 
the nuclear system, for a given electronic state. If the electronic state is also 
the ground state, this will be the ground state of the atomic system. All other 
minima are local minima, that is, metastable states that the nuclear system can 
leave by overcoming a potential barrier. 



8.3 Diatomic molecules 

Let us consider now the simple case of diatomic molecules, and in particular, 
the molecule of H2. There are 6 nuclear coordinates, 3 for the center of mass 
and 3 relative to it, but just one, the distance R between the nuclei, determines 
the effective interatomic potential V{R): the PES is in fact invariant, both 
translationally and for rotations around the axis of the molecule. Given a 



distance R, we may solve Eq.(8.4) for electrons, find the Z-th electronic energy 
level E( l >(R) and the corresponding interatomic potential V(R) = Ejj(R) + 
E^(R). Note that the nuclear repulsion energy Eu(R) is simply given by 

Eu(R) = ^ (8.11) 

where Z\ and Z<i are nuclear charges for the two nuclei. 

Let us consider the electronic ground state only for H2 molecule. At small 
R, repulsion between nuclei is dominant, V{R) becomes positive, diverging like 
q\jR for R — > 0. At large R, the ground state becomes that of two neutral 
H atoms, thus V(R) — 2Ry. At intermediate R, the curve has a minimum at 
about R = 0.74A, with V(R ) ~ V(oo) - 4.5eV. Such value - the difference 
between the potential energy of atoms at large distances and at the minimum, 
is known as cohesive or binding energy. The form of the interatomic potential is 
reminiscent of that of model potentials like Morse (in fact such potentials were 
proposed to model binding). 

What is the electronic ground state for R ~ i?o? We can get an idea by 
using the method of molecular orbitals: an approximate solution in which single- 
electron states are formed as linear combinations of atomic states centered 
about the two nuclei. Combinations with the same phase are called ligand, as 
they tend to accumulate charge in the region of space between the two nuclei. 
Combinations with opposite phase are called antiligand, since they tend to 
deplete charge between the nuclei. Starting from molecular orbitals, one can 
build the wave function as a Slater determinant. Two ligand states of opposite 
spin, built as superpositions of Is states centered on the two nuclei (a orbitals), 
yield a good approximation to the ground state. By construction, the total spin 
of the ground state is zero. 
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Molecular orbitals theory can explain qualitatively the characteristics of the 
ground (and also excited) states, for the homonuclear dimer series (i.e. formed 
by two equal nuclei). It is however no more than semi-quantitative; for better 
results, one has to resort to the variational method, typically in the framework 
of Hartree-Fock or similar approximations. Orbitals are expanded on a basis 
set of functions, often atom-centered Gaussians or atomic orbitals, and the 
Hartree-Fock or similar equations are solved on this basis set. 



8.4 Code: h2_hf_gauss 

Code h2_hf _gauss . f 9C^](or h2Jif _gauss . cQ solves the Hartree-Fock equations 
for the ground state of the H2 molecule, using a basis set of S Gaussians. The 
basic ingredients are the same as in code heliumJif _gauss, but the basis set 
is now composed of two sets of Gaussians, one centered around nucleus 1 and 
one centered around nucleus 2. As a consequence, the overlap matrix and the 
matrix element of the Fock matrix contain terms (multi- centered integrals) in 
which the nuclear potential and the two basis functions are not centered on the 
same atom. 

The code requires in input a set of Gaussians coefficients (with the usual 
format); then it solves the SCF equations at interatomic distances d m i n , d m i n + 
8d, dmin + 25d, . . ., dmax (the parameters dmin-, d m ax, Sd have to be provided 
in input). It prints on output and on file h2.out the electronic energy (not 
to be confused with the Hartree-Fock eigenvalue), nuclear repulsive energy, the 
sum of the two (usually referred to as "total energy"), all in Ry; finally, the 
difference between the total energy and the energy of the isolated atoms at 
infinite distance, in eV. 

Note that _i2J_.f _gauss . f 90 also solves the H^ case if the variable do_scf 
is set to .false.. In this case the Schrodinger equation is solved, without any 
SCF procedure. 



8.4.1 Gaussian integrals 

A complete derivation of the Gaussian integrals appearing in the solution of the 
Hartree-Fock equations can be found at pages 77-81 of the book of Thijssen. 
The following is a quick summary. 

We use a basis set of gaussian functions, in which the index now labels not 
only the coefficient of the Gaussian but also the center (one of the two nuclei 
in practice): 

6i(r) =exp(-Q i (r-R i ) 2 ) . (8.12) 
The following theorem for the product of two Gaussians: 



exp (— «j(r - Rj) 2 ) x exp (— ay(r - Rj) 2 ) = exp -(a, + a,-)(r - Rj 

(8.13) 



. ) 2 



1 http://www.fisica.uniud.it/%7Egiannozz/Corsi/MQ/Software/F90/h2_hf_gauss.f90 



http:/ /www.hsica.uniud.it/%7Egiannozz/Corsi/MQ/Software/C/h2_hf_gauss.c 
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where 



K, 



exp 



a i a j |p pi ,2 
rtj — I \ 

OLi + Oy 



R, 



OjRj + ayRj 
% + a. 



.14) 



allows to calculate the superposition integrals as follows: 

3/2 



5, 



& i (r)6 j (r)<fV 



7T 



A, 



'.)■ 



The kinetic contribution can be calculated using the Green's theorem: 



&;(r)V 2 &j(r)tfV = / V&i(r)V6j(r)d 3 r 



and finally 



OiiOL 



1^3 



ou + a, 



6 _4- a ^' 



Ri — R; 



Sij. 



.15) 



.16) 



5.17) 



The calculation of the Coulomb interaction term with a nucleus in R' ^ R is 
more complex and requires to go through Fourier transforms. At the end one 
gets the following expression: 



Vi 



bi(rh 



1 



R 



7 MT)d\ = S; 



I] 



1 

R ?1 — R 



yrerf (y/ati + ciy|Ry - R'|) . 



(8-18) 

In the case Rjj - R' = we use the limit erf(x) —> 2x/y / 7r to get the already 



known result: Vi 



9iljm 



l3 — —2-n/{ai + ctj). The bi-electronic integrals introduced in 
l be calci 



the previous chapter, Eq.(7.26), can be calculated using a similar technique: 

1 



6,,(r)6,(r) 



19) 



Sin S) 



1 



R,i — R 



7 erf I 



im I 



(ai + aj)(ai + an, 
014 + otj + ai + a Ti 



R-?' 7 R-, 



im I 



(beware indices!). 

Although symmetry is not used in the code, it can be used to reduce by a 
sizable amount the number of bi-electronic integrals gujm- They are obviously 
invariant under exchange of i, j and l,m indices. This means that if we have 
N basis set functions, the number of independent matrix elements is not A 4 
but M 2 , where M = N(N + l)/2 is the number of pairs of and (l,m) 
indices. The symmetry of the integral under exchange of r and r' leads to 
another symmetry: gnj m is invariant under exchange of the and (l,m) 

pairs. This further reduces the independent number of matrix elements by a 
factor 2, to M(M + l)/2 ~ A 4 /8. 

8.4.2 Laboratory 

• Chosen a basis set that yields a good description of isolated atoms, find the 
equilibrium distance by minimizing the (electronic plus nuclear) energy. 
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Verify how sensitive the result is with respect to the dimension of the 
basis set. Note that the "binding energy" printed on output is calculated 
assuming that the isolated H atom has an energy of -1 Ry, but you should 
verify what the actual energy of the isolated H atom is for your basis set. 

Make a plot of the ground state molecular orbital at the equilibrium 
distance along the axis of the molecule. For a better view of the binding, 
you may also try to make a two-dimensional contour plot on a plane 
containing the axis of the molecule. You need to write a matrix on a 
uniform two-dimensional N x M grid in the following format: 



x_0 


Y-0 


\psi(x. 


.0,y_0) 


x_l 


y-o 


\psi(x. 


-i,y_o) 


x_N 


y-o 


\psi(x. 


.N,y_0) 


(blank 


line) 






x_0 


y-i 


\psi(x. 


.0,y_l) 


x_l 


y-i 


\psi(x. 


-i,y_i) 


x_N 


y-i 


\psi(x. 


.N,y_l) 


(blank 


line) 






x_0 


y_M 


\psi(x. 


.0,y_M) 


x_l 


y_M 


\psi(x. 


.i,y_M) 


x_N 


y_M 


\psi(x. 


_N,y_M) 



and gnuplot commands set contour; unset surface; set view 0, 90 
followed by splot "file name" u 1:2:3 w 1 

Plot the ground state molecular orbital, together with a ligand combi- 
nation of Is states centered on the two H nuclei (obtained from codes 
for hydrogen). You should find that slightly contracted Slater orbitals, 
corresponding to Z = 1.24, yield a better fit than the Is of H. Try the 
same for the first excited state of H2 and the antiligand combination of 
Is states. 

Study the limit of superposed atoms (R — > 0) and compare with the results 
of codes hydrogen_gauss and heliumJif _gauss, with the equivalent basis 
set. The limit of isolated atoms (R — > 00) will instead yield strange 
results. Can you explain why? If not: what do you expect to be wrong 
in the Slater determinant in this limit? 

Can you estimate the vibrational frequency of H2? 
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Chapter 9 



Electrons in a periodic 
potential 

The computation of electronic states in a solid is a nontrivial problem. A great 
simplification can be achieved if the solid is a crystal, i.e. if it can be described 
by a regular, periodic, infinite arrangement of atoms: a crystal lattice. In this 
case, under suitable assumptions, it is possible to re-conduct the solution of the 
many-electron problem (really many: 0(1O 23 )!) to the much simpler problem 
of an electron under a periodic potential. Periodicity can be mathematically 
formalized in a simple and general way in any number of dimensions, but in the 
following we will consider a one-dimensional model case, that still contains the 
main effects of periodicity on electrons. 

9.1 Crystalline solids 

Let us consider an infinite periodic system of "atoms" , that will be represented 
by a potential, centered on the atomic position. This potential will in some 
way - why and how being explained in solid-state physics books - the effec- 
tive potential (or crystal potential) felt by an electron in the crystal. We will 
consider only valence electrons, i.e. those coming from outer atomic shells. 
Core electrons, i.e. those coming from inner atomic shells, are tightly bound 
to the atomic nucleus: their state is basically atomic-like and only marginally 
influenced by the presence of other atoms. We assume that the effects of core 
electrons can be safely included into the crystal potential. The pseudopotential 
approach formalizes the neglect of core electrons. 

The assumption that core electrons do not significantly contribute to the 
chemical binding and that their state does not significantly change with respect 
to the case of isolated atoms is known as frozen- core approximation. This is 
widely used for calculations in molecules as well and usually very well verified 
in practice. 

We also consider independent electrons, assuming implicitly that the crystal 
potential takes into account the Coulomb repulsion between electrons. The aim 
of such simplification is to obtain an easily solvable problem that still captures 
the essence of the physical problem. With a judicious choice of the crystal 
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potential, we can hope to obtain a set of electronic levels that can describe the 
main features of the crystal. A rigorous basis for such approach can be provided 
by Hartree-Fock or Density-Functional theory. In the end, the basic step is to 
solve to the problem of calculating the energy levels in a periodic potential. 

We haven't yet said anything about the composition and the periodicity of 
our system. Let us simplify further the problem and assume a one-dimensional 
array of atoms of the same kind, regularly spaced by a distance a. The atomic 
position of atom n will thus be given as a n = na, with n running on all integer 
numbers, positive and negative. In the jargon of solid-state physics, a is the 
lattice parameter, while the a n are the vectors of the crystal lattice. The system 
has a discrete translational invariance, that is: it is equal to itself if trans- 
lated by a or multiples of a. Called V{x) the crystal potential, formed by the 
superposition of atomic-like potentials: V(x) = ^2 n V n {x — a n ), the following 
symmetry holds: V{x + a) = V(x). Such symmetry plays a very important role 
in determining the properties of crystalline solids. Our one-dimensional space 
(the infinite line) can be decomposed into finite regions (segments) of space, of 
length a, periodically repeated. A region having such property is called unit cell, 
and the smallest possible unit cell is called primitive cell. Its definition contains 
some degree of arbitrarity: for instance, both intervals [0, a[ and ] — a/2, +o/2] 
define a valid primitive cell in our case. 

9.1.1 Periodic Boundary Conditions 

Before starting to look for a solution, we must ask ourselves how sensible it is 
to apply such idealized modelling to a real crystal. The latter is formed by a 
macroscopically large (in the order of the Avogadro number or fractions of it) 
but finite number of atoms. We might consider instead a finite system con- 
taining N atoms with N — > oo, but this is not a convenient way: translational 
symmetry is lost, due to the presence of surfaces (in our specific ID case, the 
two ends). A much more convenient and formally correct approach is to intro- 
duce periodic boundary conditions (PBC). Let us consider the system in a box 
with dimensions L = Na and let us consider solutions obeying to the condition 
yj(x) = yj(x + L), i.e. periodic solutions with period L » a. We can imagine 
our wave function that arrives at one end "re-enters" from the other side. In the 
one-dimensional case there is a simple representation of the system: our atoms 
are distributed not on a straight line but on a ring, with atom N between atom 
N — 1 and atom 1 . 

The advantage of PBC is that we can treat the system as finite (a segment 
of length L in the one-dimensional case) but macroscopically large (having N 
atoms, with N macroscopically large if a is a typical interatomic distance and 
L the typical size of a piece of crystal), still retaining the discrete translational 
invariance. Case N — > oo describes the so-called thermodynamical limit. It is to 
be noticed that a crystal with PBC has no surface. As a consequence there is no 
"inside" and "outside" the crystal: the latter is not contemplated. This is the 
price to pay for the big advantage of being able to use translational symmetry. 

In spite of PBC and of translational symmetry, the solution of the Schro- 
dinger equation for a periodic potential does not yet look like a simple problem. 
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We will need to find a number of single-particle states equal to at least half the 
number of electrons in the system, assuming that the many-body wave function 
is build as an anti-symmetrized product of single-electron states taken as spin- 
up and spin-down pairs (as in the case of He and H2). Of course the resulting 
state will have zero magnetization (5 = 0). The exact number of electrons in a 
crystal depends upon its atomic composition. Even if we assume the minimal 
case of one electron per atom, we still have N electrons and we need to calculate 
N/2 states, with N — > 00. How can we deal with such a macroscopic number 
of states? 



9.1.2 Bloch Theorem 

At this point symmetry theory comes to the rescue under the form of the Bloch 
theorem. Let us indicate with T the discrete translation operator: Tip(x) = 
ip{x + a). What is the form of the eigenvalues and eigenvectors of 7~? It can be 
easily verified (and rigorously proven) that Tip(x) = Xip(x) admits as solution 
V ; fc( a; ) = exp(i/cx)iifc(x), where A; is a real number, Uk{x) is a periodic function 
of period a: Uk(x + a) = Uk{x). This result is easily generalized to three 
dimensions, where k is a vector: the Bloch vector. States ipk are called Bloch 
states. It is easy to verify that for Bloch states the following relation hold: 

Mx + o) = MxY ka - (9.1) 

Let us classify our solution using the Bloch vector k (in our case, a one- 
dimensional vector, i.e. a number). The Bloch vector is related to the eigenvalue 
of the translation operator (we remind that H and T are commuting opera- 



tors). Eq.(9.1 ) suggests that all k differing by a multiple of 2ir/a are equivalent 
(i.e. they correspond to the same eigenvalue of T). It is thus convenient to 
restrict to the following interval of values for k: k: —ir/a < k < ir/a. Values 
of k outside such interval are brought back into the interval by a translation 
G n = 2-Kn/a. 

We must moreover verify that our Bloch states are compatible with PBC. 
It is immediate to verify that only values of k such that exp(ikL) = 1 are 
compatible with PBC, that is, k must be an integer multiple of 2ir/L. As 
a consequence, for a finite number N of atoms (i.e. for a finite dimension 
L = Na of the box), there are N admissible values of k: k n = 2irn/L, con 
n = —N/2, ...,N/2 (note that k_^/ 2 = — 7r/a is equivalent to kj^/ 2 = 7r/o). In 
the thermodynamical limit, N — > 00, these iV Bloch vectors will form a dense 
set between —ir/a and ir/a, in practice a continuum. 



9.1.3 The empty potential 

Before moving towards the solution, let us consider the case of the simplest 
potential one can think of: the non-existent potential, V(x) = 0. Our system 
will have plane waves as solutions: ipk{x) = (l/v / i)exp(ifcx), where the factor 
ensure the normalization, k may take any value, as long as it is compatible 
with PBC, that is, k = 2irn/L, with n any integer. The energy of the solution 
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with wave vector k will be purely kinetic, and thus: 



*(«) = ^f, «(*) = ?£. (9.2) 

In order to obtain the same description as for a periodic potential, we simply 
"refold" the wave vectors k into the interval — ir/a < k < ir/a, by applying 
the translations G n = 2im/a. Let us observe the energies as a function of the 



"refolded" k, Eq.(9.2): for each value of k in the interval —Tr/a < k < ir/a there 
are many (actually infinite) states with energy given by e n (k) = Ti 2 {k-\-G n ) 2 /2m. 
The corresponding Bloch states have the form 

^k,n{x) = -j^e^Uk^x), u k , n (x) = e tG " x . (9.3) 
V L 

The function u ktn (x) is by construction periodic. Notice that we have moved 
from an "extended" description, in which the vector k covers the entire space, to 
a "reduced" description in which k is limited between —ir/a and ir/a. Also for 
the space of vectors k, we can introduce a "unit cell", ] — it /a, vr/a], periodically 
repeated with period 2ir/a. Such cell is also called Brillouin Zone (BZ). It is 
immediately verified that the periodicity in /c-space is given by the so-called 
reciprocal lattice: a lattice of vectors G n such that G n ■ a rn = 2-Kp, where p is 
an integer. 

9.1.4 Solution for the crystal potential 

Let us now consider the case of a "true" , non-zero periodic potential: we can 
think at it as a sum of terms centered on our "atoms" ' : 

V{x) = ^v{x -no), (9.4) 

n 

but this is not stictly necessary. We observe first of all that the Bloch theorem 
allows the separation of the problem into independent sub-problems for each k. 



If we insert the Bloch form, Eq.(9.1), into the Schrodinger equation: 

(T + V(x))e ikx u k {x) = Ee ikx u k (x), (9.5) 
we get an equation for the periodic part u k (x): 




E 



u k (x) = (9.6) 



that has in general an infinite discrete series of solutions, orthogonal between 
them: 

fL/2 ra/2 

u* k n {x)u k m {x)dx = 5 nm N u* kn (x)u k '. m (x)dx, (9.7) 

-L/2 J-a/2 

where we have made usage of the periodicity of functions u{x) to re-conduct the 
integral on the entire crystal (from —L/2 to L/2) to an integration on the unit 
cell only (from —a/2 to a/2). In the following, however, we are not going to use 
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such equations. Notice that the solutions having different k are by construction 
orthogonal. Let us write the superposition integral between Bloch states for 
different k: 

L/2 rL/2 / 

ip* kn (x)ip k , j7n (x)dx = / e^ k -^ x u* kn (x)u k ^ m (x)dx (9.8) 

-L/2 J-L/2 

where the sum over p runs over all the N vectors of the lattice. The purely 
geometric factor multiplying the integral differs from zero only if k and k' 
coincide: 

J2e^ k '- k > = N5 k:k ,. (9.9) 
v 

we have used Kronecker's delta, not Dirac's delta, because the k form a dense 
but still finite set (there are N of them). We note that the latter orthogonality 
relation holds for whatever periodic part, u(x), the Bloch states may have. 
Nothing implies that the periodic parts of the Bloch states at different k are 
orthogonal: only those for different Bloch states at the same k are orthogonal 



(see Eq.(9.7)) 



9.1.5 Plane- wave basis set 

Let us come back to the numerical solution. We look for the solution using a 
plane- wave basis set. This is especially appropriate for problems in which the 
potential is periodic. We cannot choose any plane- wave set, though: the correct 
choice is restricted by the Bloch vector and by the periodicity of the system. 
Given the Bloch vector k, the "right" plane-wave basis set is the following: 

b n , k (x) = ie*^) 1 , G n = —n. (9.10) 
VL a 

The "right" basis must in fact have a exp(ikx) behavior, like the Bloch states 
with Bloch vector k; moreover the potential must have nonzero matrix elements 



between plane waves. For a periodic potential like the one in Eq.(9.4), matrix 
elements: 

fL/2 

(bi,k\V\hk) = T I V{x)e- lGx dx (9.11) 

-L/2 

e -i P Ga\ f a/2 V (x)e- iGx dx, (9.12) 

— / J -a/2 

where G = Gi — Gj, are non-zero only for a discrete set of values of G. In fact, 
the factor J2 P e ~ %pGa 1S zero except when Ga is a multiple of 2tt, i.e. only on 
the reciprocal lattice vectors G n defined above. One finally finds 

(bi,k\V\b j>k ) = 1 V(x)e- i ^- G ^ x dx. (9.13) 

a J -a/2 

The integral is calculated in a single unit cell and, if expressed as a sum of 
atomic terms localized in each cell, for a single term in the potential. Note that 
the factor N cancels and thus the N — > oo thermodynamic limit is well defined. 
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Fast Fourier Transform 



In the simple case that will be presented, the matrix elements of the Hamilto- 
nian, Eq.(9.13), can be analytically computed by straight integration. Another 
case in which an analytic solution is known is a crystal potential written as a 
sum of Gaussian functions: 

jV-1 

V(x) = v{x-pa), v(x) = Ae- ax \ (9.14) 

p=0 

This yields 



1 r L / 2 9 
(h,k\V\hk) = ~ / Ae~ ax e~ iGx dx (9.15) 



The integral is known (it can be calculated using the tricks and formulae given 
in previous sections, extended to complex plane): 



ax \~ iGx dx = x l-e- G2 ' ia (9.16) 



e 

(remember that in the thermodynamical limit, L — > oo). 

For a generic potential, one has to resort to numerical methods to calculate 
the integral. One advantage of the plane-wave basis set is the possibility to 
exploit the properties of Fourier Transforms (FT), in particular the Fast Fourier 
Transform (FFT) algorithm. 

Let us discretize our problem in real space, by introducing a grid of n points 
Xi = ia/n, i = 0, n — 1 in the unit cell. Note that due to periodicity, grid points 
with index i > n or i < are "refolded" into grid points in the unit cell (that 
is, V(xi +n ) = V(xi), and in particular, x n is equivalent to x$. Let us introduce 
the function fj defined as follows: 

fj = jJv{x)e- iG ^dx, G i=r^> i = 0,n-l. (9.17) 

Apart from factors, this is the usual definition of FT (and is nothing but matrix 
elements). We can exploit periodicity to show that 

/,• = - r V(x)e- iG J x dx. (9.18) 
a Jo 

Note that the integration limits can be translated at will, again due to peri- 
odicity. Let us write now such integrals as a finite sum over grid points (with 
Ax = a/n as finite integration step): 



n-l 

a 
1 



-. n— 1 

fj = - V(xm)e~ iG: > Xm Ax 



m=0 
n-l 



En 



n 

m=0 



-i n— 1 

1 7 TTt 

- V V m exp[-27r— i], VJ = V(xi). (9.19) 
n ^— n 

m=0 
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Notice that the FT is now a periodic function in the variable G, with period 
G n = 2im/a\ This shouldn't come as a surprise though: the FT of a periodic 
function is a discrete function, the FT of a discrete function is periodic. 

It is easy to verify that the potential in real space an be obtained back from 
its FT as follows: 

n-i 

V(x) = f,<' (: ''■ (9-20) 
j=0 

yielding the inverse FT in discretized form: 

71-1 

Vj = E fmexp[2ir J —i]. (9.21) 

m=0 



The two operations of Eq.(9.19) and (9.19) are called Discrete Fourier Trans 



form, or DFT (not to be confused with Density- Functional Theory!). One may 
wonder where have all the G vectors with negative values gone: after all, we 
would like to calculate fj for all j such that \Gj\ 2 < E c (for some suitably 
chosen value if E c ), not for Gj with j ranging from to n — 1. The periodicity 
of DFT in both real and reciprocal space however allows to refold the Gj on 
the "right-hand side of the box" , so to speak, to negative Gj values, by making 
a translation of 2im/a. 



9.2 Code: periodicwell 

Let us now move to the practical solution of a "true" , even if model, potential: 
the periodic potential well, known in solid-state physics since the thirties under 
the name of Kronig-Penney model: 

V(x) = y^vjx — na), v(x) = -V \x\<^, v(x) = M > ^ (9.22) 

n 

and of course a > b. Such model is exactly soluble in the limit b — > 0, Vo — >• oo, 
Vob — ^constant . 

The needed ingredients for the solution in a plane-wave basis set are almost 



all already found in Sec.M.3^ and (4.4), where we have shown the numerical 



solution on a plane-wave basis set of the problem of a single potential well. 
Code periodicwell . f9cQ (or periodicwell . cQ is in fact a trivial extension 
of code pwell. Such code in fact uses a plane- wave basis set like the one 



in Eq.(9.10), which means that it actually solves the periodic Kronig-Penney 
model for k = 0. If we increase the size of the cell until this becomes large with 
respect to the dimension of the single well, then we solve the case of the isolate 
potential well. 

The generalization to the periodic model only requires the introduction of 



the Bloch vector k. Our base is given by Eq.(9.10). In order to choose when to 



http:/ /www. fisica.uniud.it/%7Egiannozz/Corsi/MQ/Software/F90/periodicwell.f90 



http:/ /www.fisica.uniud.it/7o7Egiannozz/Corsi/MQ/Software/C/periodicwell.c 



72 



truncate it, it is convenient to consider plane waves up to a maximum (cutoff) 
kinetic energy: 

5^£<^,. (9.23) 
Bloch wave functions are expanded into plane waves: 

Mx) = J2 Cnb ^A x ) ( 9 - 24 ) 

n 

and are automatically normalized if 2~2 n \ c n\ 2 = 1- The matrix elements of the 
Hamiltonian are very simple: 

Hij = (b hk \H\b hk ) = ^ h2{k ^ m Gi)2 + - Gj), (9.25) 

where V(G) is the Fourier transform of the crystal potential. Code pwell may 
be entirely recycled and generalized to the solution for Bloch vector k. It is 
convenient to introduce a cutoff parameter E cut for the truncation of the basis 
set. This is preferable to setting a maximum number of plane waves, because 
the convergence depends only upon the modulus of k + G. The number of plane 
waves, instead, also depends upon the dimension a of the unit cell. 

Code periodicwell requires in input the well depth, Vo, the well width, 
b, the unit cell length, a. Internally, a loop over k points covers the entire BZ 
(that is, the interval [— ir/a, n/a] in this specific case), calculates E(k), writes 
the lowest E(k) values to files bands. out in an easily plottable format. 



9.2.1 Laboratory 

• Plot E(k), that goes under the name of band structure, or also dispersion. 
Note that if the potential is weak (the so-called quasi-free electrons case) , 
its main effect is to induce the appearance of intervals of forbidden energy 
(i.e.: of energy values to which no state corresponds) at the boundaries 
of the BZ. In the jargon of solid-state physics, the potential opens a gap. 
This effect can be predicted on the basis of perturbation theory. 

• Observe how E(k) varies as a function of the periodicity and of the well 
depth and width. As a rule, a band becomes wider (more dispersed, in 
the jargon of solid-state physics) for increasing superposition of the atomic 
states. 

• Plot for a few low-lying bands the Bloch states in real space (borrow and 
adapt the code from pwell). Remember that Bloch states are complex 
for a general value of k. Look in particular at the behavior of states for 
k = and k = ±ir/a (the "zone boundary"). Can you understand their 
form? 
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Chapter 10 

Pseudopotentials 



In general, the band structure of a solid will be composed both of more or less 
extended (valence) states, coming from outer atomic orbitals, and of strongly- 
localized (core) states, coming from deep atomic levels. Extended states are the 
interesting part, since they determine the (structural, transport, etc.) proper- 
ties of the solid. The idea arises naturally to get rid of core states by replacing 
the true Coulomb potential and core electrons with a pseudopotential (or effec- 
tive core potential in Quantum Chemistry parlance) : an effective potential that 
"mimics" the effects of the nucleus and the core electrons on valence electrons. 
A big advantage of the pseudopotential approach is to allow the usage of a 
plane- wave basis set in realistic calculations. 

10.1 Three-dimensional crystals 

Let us consider now a more realistic (or slightly less unrealistic) model of a 
crystal. The description of periodicity in three dimensions is a straightforward 
generalization of the one-dimensional case, although the resulting geometries 
may look awkward to an untrained eye. The lattice vectors, R n , can be written 
as a sum with integer coefficients, nf. 



of three primitive vectors, aj. There are 14 different types of lattices, known 
as Bravais lattices. The nuclei can be found at all sites + R n , where 
runs on all atoms in the unit cell (that may contain from 1 to thousands of 
atoms!). It can be demonstrated that the volume O of the unit cell is given by 
SI = ai • (a2 x as), i.e. the volume contained in the parallelepiped formed by 
the three primitive vectors. We remark that the primitive vectors are in general 
linearly independent (i.e. they do not lye on a plane) but not orthogonal. 

The crystal is assumed to be contained into a box containing a macroscopic 
number N of unit cells, with PBC imposed as follows: 



Of course, N = N± ■ iV 2 ■ N3 and the volume of the crystal is V = NQ. 



R n = mai + n 2 a 2 + n 3 a 3 



(10.1) 



V>(r + iViai + iV 2 a 2 + N 3 & 3 ) = ip(r). 



(10.2) 
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A reciprocal lattice of vectors G m such that G m • R n = 2np, with p integer, 
is introduced. It can be shown that 

G m = mibi + m 2 b 2 + m 3 b 3 (10.3) 

with nii integers and the three vectors bj given by 

27T 27T 27T fin as 

bi = — a 2 x a 3 , b 2 = — a 3 x a 2 , b 3 = — ai x a 2 (10.4) 

(note that a^ • bj = 2-nbij ) . The Bloch theorem generalizes to 

ip(r + R) = e ik ' R V>(r) (10.5) 

where the Bloch vector k is any vector obeying the PBC. Bloch vectors are 
usually taken into the three-dimensional Brillouin Zone (BZ), that is, the unit 
cell of the reciprocal lattice. 

It can be shown that there are N Bloch vectors; in the thermodynamical 
limit N — > oo, the Bloch vector becomes a continuous variable as in the one- 
dimensional case. We remark that this means that at each k-point we have 
to "accommodate" v electrons, where v is the number of electrons in the unit 
cell. For a nonmagnetic, spin-unpolarized insulator, this means v/2 filled states 
(usually called "valence" states, while empty states are called "conduction" 
states). We write the electronic states as ^k,j where k is the Bloch vector and 
i is the band index. 



10.2 Plane waves, core states, pseudopotentials 

For a given lattice, the plane wave basis set for Bloch states of vector k is 

6n,k(r) = J= e *( k + G "> r (10.6) 
V V 



where G n are reciprocal lattice vector. A finite basis set can be obtained, as 
seen in the previous section, by truncating the basis set up to some cutoff on 
the kinetic energy: 

par, 

In realistic crystals, however, E cu t must be very large in order to get a good 
description of the electronic states. The reason is the very localized character 
of the core, atomic-like orbitals, and the extended character of plane waves. 
Let us consider core states in a crystal: their orbitals will be very close to the 
corresponding states in the atoms and will exhibit the same strong oscillations. 
Moreover, these strong oscillations will be present in valence (i.e. outer) states 
as well, because of orthogonality (for this reason these strong oscillations are 
referred to as "orthogonality wiggles"). Reproducing highly localized functions 
that vary strongly in a small region of space requires a large number of delo- 
calized functions such as plane waves. 
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Let us estimate how large is this large number using Fourier analysis. In 
order to describe features which vary on a length scale S, one needs Fourier 
components up to Qmax ~ In / 5 . In a crystal, our wave vectors q — k -|- G have 
discrete values. There will be a number Np\y of plane waves approximately 
equal to the volume of the sphere of radius q ma x , divided by the volume Q, bz 
of the BZ: 

AtTO 3 87T 3 

a »s- (la8) 

The second equality follows from the definition of the reciprocal lattice. 

A simple estimate for diamond is instructive. The Is orbital of the carbon 
atom has its maximum around 0.3 a.u., so 5 ~ 0.1 a.u. is a reasonable value. 
Diamond has a " face-centered-cubic" lattice with lattice parameter ao = 6.74 
a.u. and primitive vectors: 

ai = a °{\ , \ , °) ' &2 = a ° ' &3 = a ° (°' 2~' 2~) ' ^°' 9 ^ 

The unit cell has a volume Q = ciq/4, the BZ s a volume VLbz = (27r) 3 /( a o/4)- 
Inserting the data, one finds Npy/ ~ 250, 000 plane wave, clearly too much for 
practical use. 

It is however possible to use a plane wave basis set in conjunction with 
pseudopotentials: an effective potential that "mimics" the effects of the nucleus 
and the core electrons on valence electrons. The true electronic valence orbitals 
are replaced by "pseudo-orbitals" that do not have the orthogonality wiggles 
typical of true orbitals. As a consequence, they are well described by a much 
smaller number of plane waves. 

Pseudopotentials have a long history, going back to the 30's. Born as a 
rough and approximate way to get decent band structures, they have evolved 
into a sophisticated and exceedingly useful tool for accurate and predictive 
calculations in condensed-matter physics. 



10.3 Code: cohenbergstresser 

Code cohenbergstresser . f9(j^] (or cohenbergstresser . cj^]) implements the 
calculation of the band structure in Si using the pseudopotentials published by 
M. L. Cohen and T. K. Bergstresser, Phys. Rev. 141, 789 (1966). These are 
"empirical" pseudopotentials, i.e. devised to reproduce available experimental 
data, and not derived from first principles. 



Si has the same crystal structure as Diamond: 
a face-centered cubic lattice with two atoms in 
the unit cell. In the figure, the black and red 
dots identify the two sublattices. The side of 
the cube is the lattice parameter ao- In the Di- 
amond structure, the two sublattices have the 
same composition; in the zincblende structure 
(e.g. GaAs), they have different composition. 



J http:/ /www. fisica.uniud.it/%7Egiannozz/Corsi/MQ/Software/F90/cohenbergstresser.f90 
^http:/ /www.fisica.uniud.it/%7Egiannozz/Corsi/MQ/Software/C/cohenbergstresser.c 
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d 



a 



(10.10) 



The origin of the coordinate system is arbitrary; typical choices are one of 
the atoms, or the middle point between two neighboring atoms. We use the 
latter choice because it yields inversion symmetry. The two atoms in the unit 
cell are thus at positions di = — d, d2 = +d, where 

i i r 

The lattice parameter for Si is ao = 10.26 a.u.. [[] 

Let us re-examine the matrix elements between plane waves of a potential 
V, given by a sum of spherical symmetric potentials centered around atomic 
positions: 

V(r) = EE^(l'-d M -R»l 



(10.13) 



V(r)e- iG - r dr = V Sl (G) cos(G • d), 



With some algebra, one finds: 

1 

: n Jn 

where G = Gi — Gj. The cosine term is a special case of a geometrical factor 
known as structure factor, while Vsi(G) is known as the atomic form factor: 



(10.14) 



V S i(G) = i J a V Si (r)e 



-iGr 



dr. 



(10.15) 



Cohen-Bergstresser pseudopotentials are given as atomic form factors for a few 
values of |G| , corresponding to the smallest allowed modules: G 2 = 0, 3, 4, 8, 11, .. 
in units of (27r/ao) 2 . 

The code requires on input the cutoff (in Ry) for the kinetic energy of plane 
waves, and a list of vectors k in the Brillouin Zone. Traditionally these points 
are chosen along high-symmetry lines, joining high-symmetry points shown in 
figure and listed below: 



r = 


(0,0,0), 


x = 


-(i,o,o), 

ao 


w = 


2n, 1 . 
a 2 


K = 


27T.3 3 , 
ao 4 4 


L = 


2vr^l 1 1 
ao 2 2 2 




We remark that the face-centered cubic lattice can also be described as a simple-cubic 
lattice: 

ai = a (1,0,0) , a 2 = a (0, 1, 0) , a 3 = a (0, 0, 1) 
with four atoms in the unit cell, at positions: 

'1 1 



di=a (o,-,-V d 2 = a (-,0, -V d 3 = a ( 



2' 2' 



.0 



(10.11) 

d 4 = (0,0,0). (10.12) 
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10.3.1 Laboratory 

• Verify which cutoff is needed to achieve converged results for E(k) 

• Reproduce the results in the original paper for Si, Ge, Sn. 

• Try to figure out what the charge density would be by plotting the sum 
of the four lowest wavefunctions squared at the T point. It is convenient 
to plot along the (110) plane (that is: one axis along (1,1,0), the other 
along (0,0,1) ). 

• In the zincblende lattice, the two atoms are not identical. Cohen and 
Bergstresser introduce a "symmetric" and an "antisymmetric" contribu- 
tion, corresponding respectively to a cosine and a sine times the imaginary 
unit in the structure factor: 

{bi*\V\b jtk ) = V S (G) cos(G • d) + iV a (G) sin(G • d). (10.16) 

What do you think is needed in order to extend the code to cover the case 
of Zincblende lattices? 
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Chapter 11 



Exact diagonalization of 
quantum spin models 



Systems of interacting spins are used since many years to model magnetic phe- 
nomena. Their usefulness extends well beyond the field of magnetism, since 
many different physical phenomena can be mapped, exactly or approximately, 
onto spin systems. Many models exists and many techniques can be used to 
determine their properties under various conditions. In this chapter we will deal 
with the exact solution (i.e. finding the ground state) for the Heisenberg model, 
i.e. a quantum spin model in which spin centered at lattice sites interact via 
the exchange interaction. The hyper-simplified model we are going to study is 
sufficient to give an idea of the kind of problems one encounters when trying to 
solve many-body systems without resorting to mean-field approximations (i.e. 
reducing the many-body problem to that of a single spin under an effective field 
generated by the other spins). Moreover it allows to introduce two very impor- 
tant concepts in numerical analysis: iterative diagonalization and sparseness of 
a matrix. 



11.1 The Heisenberg model 

Let us consider a set of atoms in a crystal, each atom having a magnetic mo- 
ment, typically due to localized, partially populated states such as 3d states in 
transition metals and 4/ states in rare earths. The energy of the crystal may in 
general depend upon the orientation of the magnetic moments. In many cases 
these magnetic moments tend to spontaneously orient (at sufficiently low tem- 
peratures) along a given axis, in the same direction. This phenomenon is known 
as ferromagnetism. Other kinds of ordered structures are also known, and in 
particular antiferromagnetism: two or more sublattices of atoms are formed, 
having opposite magnetization. Well before the advent of quantum mechanics, 
it was realized that these phenomena could be quite well modeled by a system 
of interacting magnetic moments. The origin of the interaction was however 
mysterious, since direct dipole-dipole interaction is way too small to justify the 

1 but not for our model: it can be demonstrated that the magnetization vanishes at T 7^ 0, 
for all 1-d models with short-range interactions only 
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observed behavior. The microscopic origin of the interaction was later found 
in the antisymmetry of the wave functions and in the constraints it imposes on 
the electronic structure (this is why it is known as exchange interaction). 

One of the phenomenological models used to study magnetism is the Heisen- 
berg model. This consists in a system of quantum spins Sj, localized at lattice 
sites i, described by a spin Hamiltonian: 

# = - E (Mij)Sx(i)S x {j) + Jy(ij)S y (i)S y (j) + J z (ij)S z {i)S z (j)) (11.1) 

<ij> 

The sum runs over all pairs of spins. 

In the following, we will restrict to the simpler case of a single isotropic 
interaction energy J between nearest neighbors only: 

H = -JJ2 S«-S(i). (11.2) 

<ij> 

The restriction to nearest-neighbor interactions only makes physical sense, since 
in most physically relevant cases the exchange interaction is short-ranged. We 
will also restrict ourselves to the case S = 1/2. 

11.2 Hilbert space in spin systems 

The ground state of a spin system can be exactly found in principle, since the 
Hilbert space is finite: it is sufficient to diagonalize the Hamiltonian over a 
suitable basis set of finite dimension. The Hilbert space of spin systems is in 
fact formed by all possible linear combinations of products: 

\H) = |<7 M (1)) ® |^(2)) ® . . . <g> MAO) (11.3) 

where ./V is the number of spins and the (tAi) labels the two possible spin states 
(<r = —1/2 or a = +1/2) for the i— th spin. The Hilbert space has dimension 
Nh = 2 N (or Nh = (2S + 1) N for spin S), thus becoming quickly intractable for 
N as small as a few tens (e.g. for N = 30, Nh ~ 1 billion). It is however possible 
to reduce the dimension of the Hilbert space by exploiting some symmetries of 
the system, or by restricting to states of given total magnetization. For a 
system of N spins, n up and N — n down, it can be easily demonstrated that 
Nh = N\/n\/{N — n)\. For 30 spins, this reduces the dimension of the Hilbert 
space to "only" 155 millions at most. The solution "reduces" (so to speak) to 
the diagonalization of the Ny t x Nh Hamiltonian matrix = (n\H\v), where 
fi and v run on all possible Nh states. 

For a small number of spins, up to 12-13, the size of the problem may still 
tractable with today's computers. For a larger number of spin, one has to resort 
to techniques exploiting the sparseness of the Hamiltonian matrix. The number 
of nonzero matrix elements is in fact much smaller than the total number of 
matrix elements. Let us re- write the spin Hamiltonian under the following form: 

H = ~\ E (S + ti)S-(j) + S-(i)S + (j) + 2S z (i)S z (j)). (11.4) 

<ij> 
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The only nonzero matrix elements for the two first terms are between states 
and \v) states such that <r M (fc) = a u (k) for all k / while for k = 

(a(i)\ m)\S + (i)S-(j)\m ® Ki)> (H.5) 

® <aO-)|5-(i)S'+(j)|a(0) ® I0O')> (H-6) 

where a(i), j3(i) mean z— th spin up and down, respectively. The term S z (i)S z (j) 
is diagonal, i.e. nonzero only for \i = v. 

Sparseness, in conjunction with symmetry, can be used to reduce the Hamil- 
tonian matrix into blocks of much smaller dimensions that can be diagonalized 
with a much reduced computational effort. 

11.3 Iterative diagonalization 

In addition to sparseness, there is another aspect that can be exploited to make 
the calculation more tractable. Typically one is interested in the ground state 
and in a few low-lying excited states, not in the entire spectrum. Calculating 
just a few eigenstates, however, is just marginally cheaper than calculating all 
of them, with conventional (LAPACK) diagonalization algorithms: an expen- 
sive tridiagonal (or equivalent) step, costing 0(iVf ) floating-point operations, 
has to be performed anyway. It is possible to take advantage of the small- 
ness of the number of desired eigenvalues, by using iterative diagonalization 
algorithms. Unlike conventional algorithms, they are based on successive re- 
finement steps of a trial solution, until the required accuracy is reached. If an 
approximate solution is known, the convergence may be very quick. Iterative 
diagonalization algorithms typically use as basic ingredients Hip, where ip is 
the trial solution. Such operations, in practice matrix-vector products, require 
0(N%) floating-point operations. Sparseness can however be exploited to speed 
up the calculation of Hip products. In some cases, the special structure of the 
matrix can also be exploited (this is the case for one-electron Hamiltonians in 
a plane- wave basis set). It is not just a problem of speed but of storage: even 
if we manage to store into memory vectors of length Nh, storage of a Nh x Nh 
matrix is impossible. 

Among the many algorithms and variants, described in many thick books, 
a time-honored one that stands for its simplicity is the Lanczos algorithm. 
Starting from \vq) = and from some initial guess \vi), we generate the following 
chain of vectors: 

K+i> = H\ Vj ) - aj\vj) - Pj\vj-i), \v j+l ) = -p—\w j+1 ), (11.7) 

where 

a-j = { Vj \H\vj), /3 j+ i = ({w j+1 \w j+ i)) 1/2 . (11.8) 

It can be shown that all vectors \vj) are orthogonal: (vi\vj) = Vz ^ j, and 
that in the basis of the \vj) vectors, the Hamiltonian has a tridiagonal form, 
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with a.j elements on the diagonal, (3j on the subdiagonal. After n steps: 
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If n = Nh, this transformation becomes exact: Ht = H, and constitutes an al- 
ternative tridiagonalization algorithm. In practice, the Lanczos recursion tends 
to be unstable and may lead to loss of orthogonality between states. If however 
we limit to a few steps, we observe that the lowest eigenvalues, and especially 
the lowest one, of matrix Ht converge very quickly to the corresponding ones 
of H. Since the diagonalization of a tridiagonal matrix is a very quick and 
easy operation, this procedure gives us a convenient numerical framework for 
finding a few low-lying states of large matrices. If moreover it is possible to 
exploit sparseness (or other properties of the matrix) to quickly calculate H\v) 
products without storing the entire matrix, the advantage over conventional 
diagonalization becomes immense. 



11.4 Code: heisenberg.exact 

Code heisenberg.exact . f 9(j^](or heisenberg.exact . finds the ground state 
energy of the 1-dimensional Heisenberg model, using Periodic Boundary Con- 
ditions: 

N 

H = -j£)S(i) • S(i + 1), S(JV + 1) = S(1). (11.10) 

i=l 

In the code, energies are in units of |J|, spins are adimensional. If J > a 
ferromagnetic ground state, with all spins oriented along the same direction, 
will be favored, while the J < case will favor an antiferromagnetic ordering. 
The sign of J is set in the code (to change it, edit the code and recompile). 

For the totally magnetized (ferromagnetic) case, the solution is trivial: there 
is just one state with all spins up (let us call it \F)), yielding Eq = (F\H\F) = 
—NJ/A. Also the case with — 1 spins up can be exactly solved. We have 
N states with N — 1 spins up, that we label as \n) = S-(n)\F). Exploiting 
translational symmetry, one introduces Bloch-like states 

1 N 

\k) = —= e lkn \n), k = 2irm/N, m = 0, N - 1. (11.11) 
v N n=l 

It can then be shown that these are the eigenvectors of H with eigenvalues 
E{k) = Eq + J(l — cos A;). Careful readers will recognize spin waves in this 
solution. A general exact solution (for other similar spin problems as well) can 
be found by means of the Bethe Ansatz, a highly nontrivial technique. 

" http://wwAv.fisica.uni ud.it/%7Egia nnozz/Co rsi/MQ/Software/F90/heisenber g_exact.f90l 
http:/ /www.fisica.uniud.it/%7Egiannozz/Corsi/MQ/Software/C/heisenberg_exact.c 
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The code requires the number N of spins and the number nup of up spins, 
computes the dimension nhil of the Hilbert space. It then proceeds to the 
labeling of the states, using a trick often employed for spin-1/2 systems: an 
integer index k, running from 1 to 2^ — 1, contains in its i—th bit the information 
(0=down, l=up) for the i—th spin. Of course this works only up to 32 spins, 
for default integers (or 64 spins for INTEGER (8)). The integer k is stored into 
array states for the states in the required Hilbert space. 

The Hamiltonian matrix is then filled (the code does not takes advantage of 
sparseness) and the number of nonzero matrix elements counted. For the S+S- 



and S-S+ terms in the Hamiltonian, only matrix elements as in 11.5 and 11.6 
respectively, are calculated. We remark that the line 

k = states(ii)+2**(j-l)-2**(i-l) 

is a quick-and-dirty way to calculate the index for the state obtained by flipping 
down spin i and flipping up spin j in state states (ii) j^] 

We then proceed to the generation of the Lanczos chain. The number nl of 
chain steps (should not exceed nhil) is prompted for and read from terminal. 
The starting vector is filled with random numbers. Note the new BLAS routines 
dnrm2 and dgemv: the former calculates the module of a vector, the latter a 
matrix- vector product and is used to calculate H\v). 

The Hamiltonian in tridiagonal form (contained in the two arrays d and e) 
is then diagonalized by the LAPACK routine dsterf , that actually finds only 
the eigenvalues. The lowest eigenvalues is then printed for increasing values of 
the dimension of the tridiagonal matrix, up to nl, so that the convergence of 
the Lanczos chain can be estimated. You can modify the code to print more 
eigenvalues. 

As a final check, the matrix is diagonalized using the conventional algorithm 
(routine dspev). Note how much slower this final step is than the rest of the 
calculation! Once you are confident that the Lanczos chain works, you can 
speed up the calculations by commenting out the exact diagonalization step. 
The limiting factor will become the size of the Hamiltonian matrix. 

11.4.1 Computer Laboratory 

• Examine the convergence of the Lanczos procedure to the exact energy 

• For the antiferromagnetic case, verify that the ground state has zero mag- 
netization. Plot the ground-state energy Eq and the first excited state E\ 
as a function of N. Try to verify if the gap E-Eq has al/JV dependence. 

• For the ferromagnetic case, verify that the ground state has all spins 
aligned. Note that the ground state will have the same energy no mat- 
ter what total magnetization you choose! This is a consequence of the 
rotational invariance of the Heisenberg Hamiltonian. Verify that the case 



with A — 1 spins up corresponds to the spin-wave solution, Eq.( 11.11). 
You will need to print all eigenvalues. 



4 A more elegant but hardly more transparent way would be to directly manipulate the 
corresponding bits. 
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Possible code extensions: 

• Modify the code in such a way that open boundary conditions (that is: the 
system is a finite chain) are used instead of periodic boundary conditions. 
You may find the following data useful to verify your results: E/N = 
-0.375 for N = 2, E/N = -1/3 for TV = 3, E/N = -0.404 per N = 4, 
E/N -»• -0.44325 per N -> oo 

• Modify the code in such a way that the entire Hamiltonian matrix is no 
longer stored. There are two possible ways to do this: 

- Calculate the Hip product "on the fly", without ever storing the 
matrix; 

— Store only nonzero matrix elements, plus an index keeping track of 
their position. 

Of course, you cannot any longer use dspev to diagonalize the Hamil- 
tonian. Note that diagonalization packages for sparse matrices, such as 
ARPACK, exist. 
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Appendix A 

Prom two-body to one-body 
problem 



Let us consider a quantum system formed by two interacting particles of mass 
mi an d m2, with no external fields. The interaction potential V(r) depends 
only upon the distance r = \r 2 — r±\ between the two particles. We do not make 
any further assumption on the form V(r). For the case of the Hydrogen atom, 
V will be of course the Coulomb potential. The Hamiltonian is 

As in the case of classical mechanics, one can make a variable change to the 
two new variables R and r: 

R = miri+m2F2 (A.2) 
mi + m 2 

r = r 2 -n (A.3) 

corresponding to the position of the center of mass and to the relative position. 
It is also convenient to introduce 

M = mi+m 2 (A. 4) 

mim 2 

m = (A.o) 

mi + m 2 

where m is known as the reduced mass. 

By introducing the new momentum operators: P = -jWr, and p = —ifbV r , 
conjugate to R and r respectively, the Hamiltonian becomes 

i.e. we have achieved separation of the variables. The center of mass behaves 
like a free particle of mass M; the solutions are plane waves. The interesting 
part is however the relative motion. The Schrodinger for the relative motion is 
the same as for a particle of mass m under a central force field V(r), having 
spherical symmetry with respect to the origin. 
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For Hydrogen (and all one-electron atoms: He + , Li 2+ , etc.,) the two parti- 
cles are a proton (or a heavier nucleus) and an electron, with a mass ratio equal 
to or larger than 1836. The reduced mass will be just a little bit smaller than 
the electron mass, m e = 0.911 x 1CT 30 Kg. 
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Appendix B 

Accidental degeneracy and 
dynamical symmetry 



The energy levels of the Coulomb potential depend only upon the main quantum 
number n. We have thus a degeneracy on the energy levels with the same n 
and different £, in addition to the one due to the 21+1 possible values of the 
quantum number m The total degeneracy (not considering spin) for a given n 
is thus 

n-l 

£(2£+l)=n 2 . (B.l) 

1=0 

The degeneracy of the energies having different i and same n is present only 
if the interaction potential is Coulombian. Although it is known as accidental 
degeneracy, it is not really "accidental". A degeneracy usually indicates the 
presence of a symmetry and thus of a conserved quantity. For instance, the 
degeneracy in m is related to spherical symmetry and to conservation of angular 
momentum. 

In classical mechanics, the equivalent of the accidental degeneracy is the 
conservation of the Runge-Lenz vector 

M = pxL _ a r 

m r 

verified for a classical Hamiltonian 

2m r' K ' ' 

This is the case of the relative motion of two bodies interacting via gravitational 
forces. In this case, the orbits are ellipses, and they are always closed: the 
orientation of the ellipses does not change with time. The Runge-Lenz vector 
lies along the major ellipses axis. The corresponding quantum vector has a 
slightly different expression: 

M = -(pxL-Lxp)-^r. (B.4) 
2m KV ^' r y ' 

It is possible to show that M is orthogonal to L, and that [M, H] = 0, i.e. it is 
a conserved quantity. 
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Appendix C 



Composition of angular 
momenta: the coupled 
representation 



Let us consider a system in which Ji and J 2 are two mutually commuting an- 
gular momentum operators. This may happen when they refer to independent 
physical systems: e.g. angular momenta of two different particles, or orbital 
and spin angular momentum of the same particle. Let us also assume that 
there are no interactions coupling them. We have four commuting observables 
that describe the system: J 2 , J\ z , J| and J2 Z - The common eigenstates are 
characterized by the set of quantum numbers: j±, mi, j'2 and m 2 . For fixed j\ 
and J2, we have thus (2ji + l)(2j2 + 1) distinct states. 

There is also another useful set of observables that describes the same sys- 
tem. We define the operator total angular momentum 



It is immediate to verify that also J satisfies to the commutation algebra of the 
angular momentum: 



since we have assumed that the commutators between components relative to 
different systems are zero. Thus [J z , J 2 ] = 0. 

We can then describe the system using the four operators J 2 , Jf , J 2 and J z . 
This is known as the coupled representation. In order to demonstrate that all 
these four operators commute, we just need to show that [J 2 , J 2 ] = [Jf , J 2 ] = 



J = Ji + J 2 



(C.l) 




= [J\x + J2x, Jly + J2y] 

= [Jlx, Jly] + [J2x, J2y] 

= iflJi z + ifl J 2z 

= ifij z 



(C.2) 
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and [JlJ z ] = [Jf , J z ] = 0: 

\J\, </ 2 ] = [J±j JxJx ± JyJy ± JzJz] 

= Jx[Jli Jx] ± \J\ i Jx]Jx + (•••) 

= Jx[J± > Jlx + <^2x] + [J\ , Jlx + J2x]Jx + (■ • •) (C-3) 

= Jx[«^l > -Ax] + [J? i «^lxRx + (• ■ ■) 

i = 

and 

[Jf, Jz] = [Jl Jiz + hz] = [Jl Ju] = (C.4) 

Let us label with j\, j'2, j and m the quantum numbers characterizing the 
eigenvalues of our operators. For fixed j\ and j will assume values between 
jmin and jmax- By definition, m = mi + mi, where m is the projection of 
Jz = Jiz + J2z- Thus, we deduce that j max = jl +32- jmm can be obtained 
by the condition that the total number of states is the same as the one for the 
original representation, i.e. 

jl+32 

J2 (2j + l) = (2j 1 + l)(2j 2 + l) (C.5) 

J — Jmin 

It can be verified that this conditions yields j m \ n = \j\ — jr 2 1 - o In summary: 
j = | ji — j2 1 , ■■ ■ ,ji +32- One can use a "vector" picture: for j = \ j\ — jr" 2 1 the two 
vectors have same direction and opposite sign, for j = j\ + 22 same direction and 
same sign, while intermediate cases correspond to angular momenta pointing 
in different directions. 

As an important example, let us consider the case j\ = 1/2 and j'2 = 1/2. 
There are four independent states, m\ = ±1/2 and 771,2 = ±1/2. Let us move 
to the coupled representation. The allowed values for j are j = and j = 1. 
For j = only m = is possible ("singlet state"). For j = 1 one has m = 0, ±1 
("triplet states"). In total, there are always four states. 

The reason to introduce the coupled representation, apparently equivalent 
to the "ordinary" one, is that in many cases there are terms in the Hamiltonian 
that couple angular momenta. Quite common for instance is the case of a 
"torque" that tends to align the two vectors: 

H = AJi -32 (C.6) 

If such term is present, J± z and J22 are no longer conserved and their eigenvalues 
are no longer good quantum numbers. In fact 

[Jiz, — AJ\ ■ J2] = —A[Ji z , Ji x J2x + JlyJ2y + J\zJ2z] (C7) 
= —A[Jlz,Jlx]J2x—A[Ji z ,Ji y ]J2y (C8) 

= -itiA(Ji y J 2x - J\xJ2y) (C.9) 

and this operator is in general not zero. 

It is however immediate to verify that J\ z + J22 is conserved, since \J-iz, — AJi- 
J 2 ] is equal to the term calculated above, but with opposite sign. 
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Appendix D 

Two-electron atoms 



Let us assume that the spin and the coordinates are separable variables (this 
is surely true if the Hamiltonian does not contain spin-dependent terms): i.e., 
one can write 

^(l,2) = $(n,r 2 )x(<7i,<7 2 ) (D.l) 

where $ is a function of coordinates r alone, \ °f the spins a alone. ^(1,2) 
is always antisymmetric because electrons are fermions. It is however clear 
that this can be achieved by an antisymmetric $ and a symmetric x, or by a 
symmetric $ and an antisymmetric x- The spin eigenfunctions of the single 
electron have two possible values that we indicate by v + and V- . We can build 
three symmetric functions for the spin: 

Xi,i = v + (a 1 )v + (a 2 ) (D.2) 

Xi,o = ^= [u+(<7i)v_(<7 2 ) + «-(<7i)t;+(o-2)] (D.3) 

Xi-i = v-{a 1 )v-{a 2 ) (D.4) 

and an antisymmetric one: 

X0,0 = ^= [v+(<7l)u_(<7 2 ) - U_(ffl)u+(<72)] (D.5) 

The symmetric functions constitute a triplet and correspond to a state of 
the two-electron system with total spin equal to 1 and three possible values: 
— 1, 0, +1, for the projection of spin along z. The antisymmetric function con- 
stitutes a singlet and corresponds to a state with total spin equal to 0. 

The value of total spin thus determines the symmetry of the spin part, and 
as of consequence, of the coordinate part of the wave function. An antisymmet- 
ric coordinate part tends to "push apart" the two electrons, since $(r, r) = 0, 
i.e. the wave function tends to zero when the two electrons move into the 
same position. The presence of electrostatic repulsion lowers the energy for 
the antisymmetric coordinate wave function with respect to the corresponding 
symmetric case, in which there is a finite probability that electrons are close 
together (no reason for <3?(r,r) to vanish). This is why excited states in He are 
labeled as ortho-helium (triplet state with symmetric spin part and antisym- 
metric coordinate part) and •para-helium (singlet state with antisymmetric spin 
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part and symmetric coordinate part), with the former lower in energy than the 
latter for the same single electron levels. The ground state of He turns out to 
be a singlet, so the coordinate wave function must be symmetric. 



D.l Perturbative Treatment for Helium atom 

The Helium atom is characterized by a Hamiltonian operator 

2m e r\ 2m e T2 rj.2 

where ryi = ~ r i\ is the distance between the two electrons. The last term 
corresponds to the Coulomb repulsion between the two electrons and makes the 
problem non separable. 

As a first approximation, let us consider the interaction between electrons: 

V=^- (D.7) 
as a perturbation to the problem described by 

ffo _|v|_ M _|v|_ M 

2m e r\ 2m e T2 

which is easy to solve since it is separable into two independent problems of 
a single electron under a central Coulomb field, i.e. a Hydrogen-like problem 
with nucleus charge Z = 2. The ground state for this system is given by the 



wave function described in Eq.(2.29) (Is orbital): 



73/2 

b°(Ti) = —=e~ Zr * (D.9) 

V7T 



in a.u.. We note that we can assign to both electrons the same wave function, as 
long as their spin is opposite. The total unperturbed wave function (coordinate 
part) is simply the product 

7 3 

,/,°( ri ,r 2 ) = _ e - z ( ri+r2 ) (D.10) 

7T 

which is a symmetric function (antisymmetry being provided by the spin part). 
The energy of the corresponding ground state is the sum of the energies of the 
two Hydrogen-like atoms: 

E = -2Z 2 Ry = -8Ry (D.ll) 

since Z = 2. The electron repulsion will necessarily raise this energy, i.e. make 
it less negative. In first-order perturbation theory, 

E-E Q = (Vol^lV'o) (D.12) 

= % I A e -2Z(r 1+ r 2 ) d 3 ri(i 3 r2 (m3) 

= \zRy. (D.14) 
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For Z = 2 the correction is equal to 2.5 Ry and yields E = —8 + 2.5 = —5.5 Ry. 
The experimental value is — 5.8074Ry. The perturbative approximation is not 
accurate but provides a reasonable estimate of the correction, even if the "per- 
turbation", i.e. the Coulomb repulsion between electrons, is of the same order 
of magnitude of all other interactions. Moreover, he ground state assumed in 
perturbation theory is usually qualitatively correct: the exact wave function for 
He will be close to a product of two Is functions. 



D.2 Variational treatment for Helium atom 

The Helium atom provides a simple example of application of the variational 
method. The independent-electron solution, Eq.(D.lO), is missing the phe- 
nomenon of screening: each electron will "feel" a nucleus with partially screened 
charge, due to the presence of the other electron. In order to account for this 
phenomenon, we may take as our trial wave function an expression like the 
one of Eq.(D.lO), with the true charge of the nucleus Z replaced by an "ef- 
fective charge" Z e , presumably smaller than Z. Let us find the optimal Z e 
variationally, i.e. by minimizing the energy. We assume 



V>(ri,r 2 ; Z e ) 



and we re- write the Hamiltonian as: 



Z e c -z e (n+n) 

7T 



(D.15) 



H 



h 2 vj 



1 ft 2 v| 



2m e ri 
We now calculate 



2m P 



T 2 



+ 



{Z-Z e )ql {Z-Z e )ql 



r-2 



+ 



r\2 
(D.16) 



E{Z e )= J iP*(r u r 2 ;Z e )H^(r 1 ,r 2 ;Z e )d 3 r 1 d 3 r 2 (D.17) 

The contribution to the energy due to the first square bracket in Eq.(D.16) is 
— 2Zg a.u.: this is in fact a Hydrogen-like problem for a nucleus with charge Z e , 
for two non-interacting electrons. By expanding the remaining integrals and 
using symmetry we find 



E(Z e ) = - 
(in a.u.) with 



-2Zi 



ri 



|Vf — d 3 rid 3 r 2 
ri2 



« -2Z e (n+r 2 ) 



7T^ 



Integrals can be easily calculated and the result is 



E(Z e 



-2Z 2 e - 4(Z - Z e )Z e + 2-Z e 



n 27 

2Zl - —Z, 



(D.18) 



(D.19) 



(D.20) 



where we explicitly set Z = 2. Minimization of E{Z e ) with respect to Z e 
immediately leads to 



27 
16 



1.6875 



(D.21) 
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and the corresponding energy is 

729 

E = ~— = -5.695 Ry (D.22) 

This result is definitely better that the perturbative result E = — 5.50Ry, 
even if there is still a non-negligible distance with the experimental result 
E = -5.8074 Ry. 

It is possible to improve the variational result by extending the set of trial 



wave functions. Sect. (7.1) shows how to produce the best single-electron func- 
tions using the Hartree-Fock method. Even better results can be obtained 
using trial wave functions that are more complex than a simple product of 
single-electron functions. For instance, let us consider trial wave functions like 

V>(ri,r 2 ) = [/(ra)<?(r 2 ) + 5 (ri)/(r 2 )] , (D.23) 

where the two single-electron functions, / and g, are Hydrogen- like wave func- 



tion as in Eq.(D.9) with different values of Z, that we label Zf and Z g . By 
minimizing with respect to the two parameters Zf and Z g , one finds Zf = 2.183, 
Z g = 1.188, and an energy E = —5.751 Ry, much closer to the experimental 
result than for a single effective Z. Note that the two functions are far from 
being similar! 

D.3 "Exact" treatment for Helium atom 

Let us made no explicit assumption on the form of the ground-state wave func- 
tion of He. We assume however that the total spin is zero and thus the coordi- 
nate part of the wave function is symmetric. The wave function is expanded over 
a suitable basis set, in this case a symmetrized product of two single-electron 
gaussians. The lower-energy wave function is found by diagonalization. Such 
approach is of course possible only for a very small number of electrons. 

Code helium_gauss . f 9cQ (or helium.gauss . looks for the ground state 
of the He atom, using an expansion into Gaussian functions, already introduced 
in the code hydrogen-gauss. We assume that the solution is the product of a 
symmetric coordinate part and of an antisymmetric spin part, with total spin 
S = 0. The coordinate part is expanded into a basis of symmetrized products 
of gaussians, B^: 

^(ra,r 2 ) = 2c fc S fc (n ) r 2 ). (D.24) 

k 



If the b{ functions are S-like gaussians as in Eq.(5.22), we have: 



B k (r 1 ,r 2 ) = -j= (&i(fc)(ri)6j(jfc)(r 2 ) + 6i(jfc)(r 2 )6 3 -(fc)(ri)) (D.25) 

where k is an index running over n(n+l)/2 pairs i(k),j(k) of gaussian functions. 
The overlap matrix Skk' may be written in terms of the Sij overlap matrices, 



Eq.(5.25), of the hydrogen-like case: 

Skk' = {B k \Bk') = {Su'Sjj' + Sij'Sji') . (D.26) 



J http:/ /www. fisica.uniud.it/%7Egiannozz/Corsi/MQ/Software/F90/helium_gauss.f90 



http:/ /www.fisica.uniud.it/%7Egiannozz/Corsi/MQ/Software/C/helium_gauss.c 
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The matrix elements, Hkk', of the Hamiltonian: 
Hkk' = (Bk\H\Bk>) , H - 



U 2 V 2 Zq\ h 2 V 2 Zql 



2m P 



n 



2m P 



+ 



r 2 ri2 



(D.27) 



+ HYa, obtained for the 

Li L J 



can be written using matrix elements Hy 
hydrogen-like case with Z = 2, Eq.(5.26) and ( 5.27| ): 

Hkk' = {Hu'Sjji + HijiSji/ + Sii'Hjjt + Hij'Sji/) + (Bk\V ee \Bk') , (D.28) 

and the matrix element of the Coulomb electron-electron interaction V ee : 



(Bk\ V ee \Bk' 



q 2 

h(k) ( r i)bj( k ) O2) -^-k(k') ( r i)^'(fc') (r 2 )d 3 rid 3 r 2 (D.29) 

f Q 2 s s 

+ / (ri)6j(fc) C r 2) — ^(fc') (ri)^(fc')( r 2)^ nd r 2 . 



+ 



a/3(a + /3) 1 /2 a '/3'(a' + /3') 1/2 ' 



(D.30) 



These matrix elements can be written, using Eq.(|7.33|), as 

(Bk\V ee \Bk>) 
where 

« = Oj(fc) + «i(fc')> /5 = %(fc) + aj(fe')' a ' = a «(fc) + a j(fc')> = aj(fc) + «i(fc')- 

(D.31) 

In an analogous way one can calculate the matrix elements between sym- 
metrized products of gaussians formed with P-type gaussian functions (those 
defined in Eq. 5.23). The combination of P-type gaussians with L = has the 
form: 



J3fc(ri,r 2 ) = ^( r i - r 2) (hk)( r i) b j(k)( r 2) + 6 i(fc) (r 2 )6 j(fc) (ri) 



(D.32) 



It is immediately verified that the product of a S-type and a P-type gaussian 
yields an odd function that does not contribute to the ground state. 

In the case with S-type gaussians only, the code writes to file "gs-wfc.out" 
the function: 

P(n,r 2 ) = (47rrir 2 ) 2 |^(n,r 2 )| 2 , (D.33) 

where P(ri, r 2 )dridr 2 is the joint probability to find an electron between r\ and 
T\ + dr\ , and an electron between r 2 and r 2 + dr 2 . The probability to find an 
electron between r and r + dr is given by p(r)dr, with 

p(r) = Airr 2 J \ip(r,r 2 )\ 2 A-Krldr 2 = J P(r,r 2 )dr 2 . (D.34) 

It is easy to verify that for a wave function composed by a product of two 
identical functions, like the one in (D.10), the joint probability is the product 
of single-electron probabilities: P(r\,r 2 ) = p{r\)p{r 2 ). This is not true in 
general for the exact wave function. 
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D.3.1 Laboratory 

• observe the effect of the number of basis functions, to the choice of coef- 
ficients A of the gaussians, to the inclusion of P-type gaussians 

• compare the obtained energy with the one obtained by other methods: 



perturbation theory with hydrogen-like wave functions, (Sec D.l), varia- 



tional theory with effective Z (Sec|D.3), exact result (-5.8074 Ry). 



Make a plot of the probability P{r\,T2) and of the difference P[r\,r2) 
p{ r i)p{ r 2) > using for instance gnuplot and the following commands: 



set view 0, 90 
unset surface 
set contour 

set cntrparam levels auto 10 

splot [0:4] [0:4] "gs-wfc.out" u 1:2:3 w 1 



Note that the probability P(r\,r2) (column 3 in "splot") is not exactly 
equal to the product p(r\)p{r2) (column 4; column 5 is the difference 
between the two). 
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Appendix E 

Useful algorithms 



E.l Search of the zeros of a function 

Short notes on the problem of numerical determination of the zeros of a function, 
i.e. the values of x that solve f{x) = for a given function. 
There are two different types of zeros: 

• odd - when f(x) changes sign at xo 

• even - when f(x) does not change sign at xo 

The latter case is more difficult from a numerical point of view. Imagine to 
add a small numerical noise to f(x): the zero disappears, or else it splits into 
two odd zeros. It is usually convenient to resort to methods to localize extrema 
(minima or maxima) rather than trying to pinpoint the zero. In the following 
we will deal with odd zeros only. 

E.l.l Bisection method 

We first have to find an interval [a, b] that includes a single zero, in such a way 
that f(a)f(b) < 0. The bisection algorithm halves the interval at each iteration, 
by refining the estimate of xq: 

1. c= (a + b)/2 

2. if /(a)/(c) < we re-define b = c; else if f(b)f(c) < we re-define a = a 

3. In this way we get a new interval [a, b] of half the width at the previous 
step, and we repeat the procedure. 

Convergence is guaranteed (it is impossible to "miss" the zero), and the loga- 
rithm of the uncertainty on the solution linearly decreases with the number of 
iterations. Some difficulty may arise in the choice of the stopping criterion. For 
instance: 

• Absolute error: \a — b\ < e. This may run into trouble if xo is very large: 
rounding errors \a — b\ might be larger then e. 
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• Relative error: \a — b\ < ea. This may run into trouble close to x = 0. 

• If the slope of f{x) close to zero is very small, there might be a finite 
interval in which f{x) is indistinguishable from zero in machine represen- 
tation. 

E.1.2 Newton-Raphson method 

The function is linearly approximated at each iteration in order to obtain a 
better estimate of the zero. Let us assume that we know both f(x) and f'(x). 
Then, close to x, 

f(x + 5)~f(x) + f'(x)5 (E.l) 

and thus, to first order, 

'-m (E - 2 > 

would yield f(x + 5) = 0. We proceed by iterating in this way. It is possible to 
show that the rate of convergence is quadratic, i.e. the number of significant 
figures approximately doubles at each iteration (while with bisection method it 
grows linearly). 

The problem of this method is that convergence is not guaranteed, notably 
when f'{x) varies a lot close to the zero. Moreover, the method assumes that 
f'{x) can be directly calculated for any given x. In cases in which this snot 
true and the derivative must be calculated via finite differences, it is preferable 
to use the secant method described below. 

E.1.3 Secant method 

This is based on a linear expansion of f(x) between two successive points, x n e 
x n+ i, of an iteration: 

f(x) = /(x n _i) + X ~_ Xn ~ X [/(*„) - f(x n -i)] (E.3) 

that provides as an estimate for the zero 

X n+1 = X n - f[X n )— ; r r (E.4) 

f{Xn) ~ f(Xn-l) 

The procedure consists in iterating the above step. It is not necessary that 
the zero is contained inside the considered interval. This however may lead to 
non-convergence in some pathological cases. In regular cases, the speed of con- 
vergence is much better than for the bisection method, although slightly slower 
that the Newton-Raphson method (which however requires the knowledge of 
the derivative). 

In most difficult cases, it is convenient to split the search into two step: first 
bisection, with a sure and safe bracketing of the zero; then the secant method 
that quickly stabilizes the searched value up to the required precision. 
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